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The  study  reported  in  this  dissertation  investigates  the  effect  of  surface  tension  on 
film  flows  on  rollers  in  the  presence  of  stationary  objects  due  to  its  relevance  in 
hydrodynamic  and  elastohydrodynamic  lubrication.  An  analysis  is  introduced  to  predict 
the  film  layer  thickness  in  the  free  surface  as  an  extension  of  the  Landau  and  Levich 
theory  ^lnd  is  utilized  in  the  study.  Parameters  of  liquid  entrainment  are  considered  to 
solve  the  numerical  problem  and  to  examine  the  effect  of  the  boundary  conditions 
characterizing  the  formation  of  the  liquid  film.  It  is  shown  in  the  analysis  that  surface 
tension,  gravity  and  inertia  effects  have  considerable  influence  in  determining  the  film 
thickness  and  its  associated  boundary  conditions. 

A mathematical  model  is  also  developed  for  parametric  investigation  of  the 
influence  of  the  rheological  properties  of  the  lubricant  in  the  thermohydrodynamic  film 
conditions  which  occur  in  slider  and  journal  bearings.  It  is  shown  in  the  analysis  that 
yield  in  shear  have  considerable  influence  on  the  flow  velocity  distribution  and  the 
pressure  developed  in  thin  film  lubrication.  A parametric  investigation  is  performed  for 
predicting  the  thickness  of  the  shear  zone  where  heat  generation  occurs  in  fluid  films  with 
fixed  geometry  between  the  stationary  and  moving  surfaces.  The  investigation  shows 
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that  the  shear  zone  thickness  is  significantly  smaller  than  the  fluid  film  thickness  in  the 
lubrication  zone.  The  results  of  this  procedure  also  show  the  same  trend  as  the 
experimentally  determined  square  root  relationship  between  the  pressure  and  sliding 
velocity  in  journal  and  slider  bearings. 
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CHAPTER  1 
INTRODUCTION 


Many  applications  in  engineering  practice  involve  the  entrainment  of  liquid  films 
on  solid  objects.  To  identify  some  of  these  systems,  these  applications  may  be  classified 
as  problems  in  withdrawal,  removal,  and  drainage  [1,2]. 

In  withdrawal  the  motion  of  an  object  from  a liquid  container  causes  entrainment 
of  the  wetting  liquid.  Withdrawal  is  limited  to  the  case  where  a part  of  the  solid  support 
remains  in  contact  with  the  liquid  container.  In  removal,  an  object  of  finite  length  is 
entirely  removed  from  a liquid  container  so  that  the  entrained  film  is  separated  from  the 
container.  In  drainage,  when  the  solid  support  is  stationary  the  liquid  flows  downward  by 
gravity.  These  three  problems  are  closely  related. 

In  these  problems,  the  mass  flow  rate  or  flux  of  the  liquid  is  unknown.  Thus  the 
class  of  problems  considered  here  is  quite  different  from  that  of  liquid  film  situations 
where  the  flux  is  known,  such  as  in  flow  down  wetted  wall.  The  flux  is  closely  related  to 
the  entrainment,  the  film  distribution,  and  the  total  amount  of  liquid  removed. 

To  the  engineer,  the  most  important  parameter  is  probably  the  flux  of  the  liquid. 
Other  parameters  of  interest  are  the  thickness  of  the  film  and  the  velocity  distribution  in 
the  film.  The  engineer  would  like  to  predict  of  these  quantities  as  functions  of  position, 
time,  and  the  physical  properties  of  the  liquid.  Knowledge  of  the  dynamics  of  these  films 
is  useful  in  many  practical  fields  such  as  coating,  cleaning,  draining,  and  lubrication.  The 
problems  as  indicated  above  are  liquid  flows  and  show  important  roles  of  fluid  dynamics. 

The  phenomenon  of  fluid  film  has  been  of  considerable  interest  in  areas  such  as 
polymer  engineering.  In  the  tribological  problems  it  also  has  been  an  active  subject  for 
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research.  There  is  a wide  range  of  applications  where  these  problems  arise  as  indicated  in 
Table  1. 


Table  1 : General  Application 


Application 

Specific  application 

Classification 

Coating 

film  and  wire 

withdrawal 

photographic  film 

Cleaning 

electroplating 

removal 

rinsing 

Draining 

drainage 

drainage 

Lubrication 

gears 

withdrawal 

bearings 

■disk 

1 . 1 Outline  of  the  Problem 


The  objective  of  this  study  is  to  develop  a mathematical  model  for  thin  film 
lubrication  problems  which  occur  in  practice.  The  first  consideration  is  to  model  the 
hydrodynamic  phenomenon  between  a rotating  disc  and  a fixed  surface  in  Newtonian 
fluids.  The  second  consideration  is  the  thermohydrodynamic  phenomenon  of  lubricating 
films  between  rotating  (or  sliding)  solids  and  any  adjacent  constraint.  This  also  includes 
the  fluid  film  lubrication  that  is  associated  with  elements  in  relative  motion. 

The  first  part  deals  with  the  flow  of  a viscous  fluid  on  the  surface  of  a partially 
submerged,  rotating  bearing  (disk)  which  is  one  of  the  basic  coating  (free  surface) 
problem  as  shown  in  part  A of  Figure  1-1 -(a)  and  (b).  The  main  features  of  this  type  of 
flow  for  a Newtonian  fluid  are  described  by  the  Landau-Levich  theory  [3].  The  problem 
is  treated  using  a nonlinear  approach  without  the  inertia  effects.  In  the  free  surface  flow 
the  formation  of  the  liquid  films  is  characterized  by  a meniscus  thickness  decreasing  with 
distance  above  the  container  and  eventually  approaching  a constant  thickness.  The 
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resultant  constant  thickness  is  a function  of  the  velocity  of  withdrawal,  gravity,  and  the 
physical  properties  of  liquid  (density,  viscosity,  and  surface  tension).  The  meniscus 
region  is  characterized  at  the  interface  by  a stagnation  point  [4].  The  film  maybe  divided 
about  this  point  into  thin  and  thick  menisci.  Here  such  parameters  of  liquid  entrainment 
occurring  from  a container  will  be  considered.  Thus,  the  main  point  for  solving  this 
problem  is  the  boundary  conditions  characterized  by  the  formation  of  the  liquid  film  [3, 
5,  6,  7].  This  can  make  two  different  flow  regions  match  at  a common  boundary. 
Without  the  use  of  asymptotic  techniques  or  numerical  methods  such  as  the  finite  element 
method  it  will  be  possible  to  solve  the  full  Navier-Stokes  equations  based  on  this  key 
point.  The  model  gives  a relationship  between  the  dimensionless  film  thickness,  T,  and 
the  Capillary  number,  Ca.  To  develop  valid  viscous  free  surface  flows  at  relatively  low 
Reynolds  numbers,  a boundary-layer  type  of  approximation  is  introduced  [8-15]. 

When  a journal  or  a disk  rotates  about  its  axis  while  partially  submerged  in  a 
liquid,  a layer  of  the  liquid  is  entrained  and  withdrawn  from  the  liquid  container.  In  the 
absence  of  a constraining  boundary  the  film  thickness  is  found  to  exhibit  a strong 
dependence  on  disk  speed,  fluid  viscosity  and  surface  tension.  In  the  free  surface  flows 
with  stationary  constraint  the  formation  of  these  liquid  films  is  characteristically  the 
hydrodynamic  lubrication  problem  as  shown  in  part  B of  Figure  1-1 -(b).  For  lightly 
loaded  situations  wide  ranging  numerical  approaches  are  carried  out  to  predict  the  film 
rupture  location  and  pressure  distributions  using  the  cavitation  boundary  conditions  of 
Coyne  and  Elrod  and  two  additional  Swift-Stieber  (Reynolds)  and  separation  boundary 
models  [16-21].  For  over  a century  the  Reynolds  equation  has  been  used  to  analyze  this 
problem.  When  the  lubrication  film  thickness  profile  is  known  and  under  idealized 
conditions  of  isothermal,  highly  viscous  Newtonian  flows  the  basic  theory  of  lubrication 
is  generally  expressed  by  the  Reynolds  equation  [22].  Classical  lubrication  mechanics 
constitutes  the  solutions  to  Reynolds  equation  with  various  boundary  conditions 
simulating  various  system  configurations.  However,  much  of  the  modem  lubrication 
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literature  deals  with  departures  which  are  due  to  thermal  effects,  inertial  effects,  and 
undetermined  film  thickness  profiles.  Consequently,  this  problem  will  be  dealt  with  as  an 
extension  of  the  first  procedure  with  the  boundary  conditions  and  parameters  of  liquid 
entrainment  occurring  from  a layer  of  the  liquid  on  the  rotating  system  entrained  and 
withdrawn  from  the  container.  While  the  film  thickness  in  the  absence  of  a constraining 
boundary  seems  to  have  a strong  dependence  on  the  rotating  speeds,  as  well  as  fluid 
viscosity  and  interfacial  tension,  the  effects  of  an  adjacent  constraining  boundary  requires 
a mathematical  model  based  on  lubrication  theory.  This  study  will  be  concerned  with  the 
hydrodynamic  behavior  for  the  critical  conditions  at  the  onset  of  film  formation 

The  last  procedure  deals  with  the  development  of  a thermohydrodynamic 
computational  procedure  for  evaluating  the  pressure,  temperature  and  velocity 
distributions  in  thin  fluid  films  with  fixed  geometry  between  the  stationary  and  moving 
bearing  (disk)  surfaces  as  shown  in  part  C of  Figure  1-1 -(b).  When  temperature  rise 
across  the  system  such  as  bearing  is  small,  the  classical  isothermal  theory  with  an 
effective  viscosity  may  be  applied.  However,  in  practical  cases,  thermal  effects  may  be 
significant,  thus,  the  classical  theory  is  no  longer  adequate  [23].  The 
thermohydrodynamic  (THD)  theory  of  lubrication  is  applied  to  evaluate  local  viscosity. 
Further,  in  most  industrial  applications  where  the  rotating  speed  is  high,  the  systems  are 
often  operated  in  a special  regime.  In  practical  applications,  when  the  lubricant  is 
subjected  to  high  shear  force  which  act  for  only  a relatively  short  period  of  time,  the 
assumption  that  the  entire  film  behaves  as  a Newtonian  fluid  may  not  be  valid.  A 
primary  hypothesis  of  this  study  is  that  in  fluid  film  lubrication  there  exists  a small  shear 
zone  at  the  middle  thickness  of  the  film  where  the  lubricant  undergoes  a high  shear  rate 
[24-27].  This  shear  zone  thickness  inside  the  lubricant  film  is  controlled  by  the  surface 
tension  at  the  solid  boundaries  and  the  rheological  behavior  of  the  lubricant  under 
different  pressure,  temperature,  and  speed  conditions.  The  velocity  variations  and  the 
heat  generation  are  assumed  to  occur  intensively  in  a central  zone  having  rapidly 
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changing  high  shear  rate  with  the  same  length  (or  length  of  the  lubrication  zone  obtained 
from  the  second  procedure)  and  width  as  the  bearing  (disk)  but  with  a significantly 
smaller  thickness  than  the  fluid  film  thickness.  A transient  thermohydrodynamic  model 
with  a rectangular  computational  domain  is  utilized  which  generalizes  the  analysis  for 
any  given  film  geometry.  The  developed  procedure  also  gives  an  analytical  basis  for 
explaining  the  “Fogg  effecf’  (thermal  wedge)  where  significant  pressures  can  be 
generated  in  bearings  as  a result  of  the  thermal  expansion  of  the  film  in  the  direction  of 
the  thickness  [28]. 

A numerical  solution  to  the  thermohydrodynamic  lubrication  with  consideration 
of  both  the  thermal  and  non-Newtonian  effects  have  been  obtained.  The  computational 
technique  requires  the  simultaneous  solution  of  the  Newtonian  hydrodynamic,  the  load, 
the  viscosity,  and  the  energy  equations.  Most  information  for  computation  of  this  part 
come  from  the  two  previous  procedures.  Some  parameters  are  involved  in  determining 
the  rotating  system  performance,  such  as  bearing,  disk.  The  analysis  on  a case-by-case 
basis  is  performed  for  the  design  of  such  as  journal  bearings,  finite  slider  bearings  and 
rotating  cylinder  (disc)  - plate  systems. 

1.2  Objective 

The  main  objective  of  the  present  study  is  to  investigate  the  role  of  surface  tension 
in  forming  adhered  films  for  tribology  and  attempt  to  theoretically  predict  the  dimension 
of  the  shear  zone  in  thermohydrodynamic  lubrication  by  using  a rheological  model  for  the 
lubricant  which  is  expected  to  behave  rheologically  under  pressure.  The  present  study 
searches  for  functions  that  approximate  the  cavitation  point  and  the  distribution  of  film 
pressure,  temperature,  and  velocity. 
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Figure  1.1  Geometry  of  simple  rotating  systems: 

(a)  Freely  confined  withdrawal 

(b)  Partially  confined  withdrawal 


CHAPTER  2 

LUBRICATION  FLUID  MECHANICS 


2.1  Introduction 


The  subject  of  hydrodynamics  is  concerned  with  the  motion  of  fluids.  There  are 
various  types  of  fluids.  Essentially  the  study  of  hydrodynamics  is  the  study  of  rates  of 
flow  of  fluid  in  conjunction  with  momentum  of  the  fluid  and  its  action  under  the  forces 
placed  upon  it. 

The  study  of  the  motion  of  fluids  is  very  important  to  lubrication,  hence 
hydrodynamics  is  a basic  element  of  lubrication  theory.  However,  hydrodynamics  alone 
is  insufficient  to  describe  the  full  characteristics  of  lubrication  because  the  lubricant 
properties  are  temperature  dependent  [29-32].  For  example,  the  viscosity  of  a lubricant  is 
a dissipative  factor  turning  mechanical  energy  into  heat.  Thus  many  applications  in 
lubrication  must  consider  the  thermal  effects  which  are  connected  with  the  viscosity  in 
the  hydrodynamics  of  the  system.  Another  important  element  of  lubrication  theory  is  the 
relationship  between  hydrodynamics  and  the  flow  passage  [29].  The  geometry  of  the 
flow  passage  can  change  with  pressure.  This  greatly  couples  hydrodynamics  with 
deformation  mechanics  in  some  cases  such  as  the  elastohydrodynamic  lubrication.  This 
phenomenon  becomes  a very  important  design  factor  in  systems  where  the  pressure  in  the 
lubricating  film  is  high  enough  to  deform  the  solids. 

The  basic  equations  for  a viscous  fluid  are  the  Navier-Stokes  equations.  They  are 
derived  by  applying  to  a fluid  the  principle  of  conservation  of  momentum.  Although  the 
general  Navier-Stokes  equations  are  rather  complicated,  approximations  can  be  made  for 
the  study  of  thin  film  of  fluids.  The  thin  film  approximation  is  the  classic  approximation 
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for  boundary  layer  theory.  However,  in  addition  to  the  thin  film  approximation,  in  most 
situations  the  motion  of  fluids  in  thin  film  is  so  viscous  that  the  momentum  convected  by 
the  fluids  is  insignificant.  Consequently,  inertial  terms  can  be  discarded  in  the  Navier- 
Stokes  equations  which  results  in  the  Reynolds  equation,  which  is  the  nucleus  of  the 
hydrodynamic  theory  of  lubrication. 

In  steady  state  and  two  dimensional  flow,  the  Navier-Stokes  equations  reduces  to 
a pair  of  elliptic  partial  differential  equations  as  follows. 
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2.2  Lubrication  Approximation  (Boundary  Laver  Theory) 


Considering  a thin  liquid  film  on  a flat  substrate,  as  shown  in  Figure  2.1,  the  flow 
is  two-dimensional  with  B,  a characteristic  length  in  the  x direction  and  h,  a characteristic 
length  in  the  z direction  [33].  Now  using  the  following  dimensionless  variables  to  scale 
the  Equations  (2.1),  (2.2),  and  (2.3), 

u=Uou‘  , x = Bx  , z = hz 
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Then  Equation  (2.3)  becomes 
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Equation  (2.1)  becomes 


9 


f 2 • 2 

Uo  • Uq  • OU 

u — r + ~T"W 


A 


V 


B 3c  B ^ 


J 


MUo  . Uo^ 

3c*  ^ 


1 * 
d u 


dz 


*2 


+ 


h 

kBj 


dx 


•2 


+ tL^s 


(2.6) 


Then,  introducing  the  assumption  for  thin  film,  h/B«l,  and  the  Reynolds  number. 


Re=puoB//d,  the  above  Equation  (2.6)  becomes 
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Now  applying  similarly  Equation  (2.2)  becomes 
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In  general,  for  many  lubrication  problems  Re(h/Bp«l,  and  consequently  the  inertia 
force  terms  of  the  left  hand  side  of  Equation  (2.7)  can  be  neglected.  However,  when 
Reynolds  number  is  large  Re(h/Bp«l  must  be  retained  in  Equation  (2.7).  However, 
Re(h/B)"^«l  is  negligible,  so  that  the  left  hand  side  of  Equation  (2.8)  is  essentially  zero. 
Since  (h/B)^«l,  Equation  (2.8)  simply  becomes 


dz 


= 0 


(2.9) 


Applying  this  to  the  three  dimensions  case,  the  dimensional  form  of  the  equations 


for  this  problem  becomes  as  follows. 


3c 


d^u 


\ 


+pg 


(2.10) 


dy 


= /^ 


(2.11) 


dz 


0 


(2.12) 


3i  dv  dw 

dx  ^ dz  ^ 


(2.13) 
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Neglecting  the  inertia  and  gravity  terms  in  Equation  (2.10)  throughout  Equation  (2.13) 
are  the  basic  equations  for  the  lubrication  approximation.  Because  of  their  very  simple 
structure  they  can  be  solved  directly  and  expressed  in  alternative,  more  useful  forms.  The 
following  form  is  for  the  x-component  of  velocity,  Equation  (2.10),  for  any  point  in  the 
lubricant  film 

» = + + (2.14) 

2 /u  ax 

v = \—^z^  + D\Z  + D2  (2.15) 

2 ju  ay 

where  C,,  Cj,  D,,  and  Dj  are  constants.  The  physical  significance  of  equation 
(2.14)  is  that  it  consists  of  two  flow  components  of  velocity,  which  represent  the 
hydraulic  flow  and  the  shear  flow,  within  the  film. 

Equation  (2.14)  is  simply  the  equation  describing  the  velocity  distribution 
between  two  surfaces  with  an  imposed  pressure  gradient.  This  is  a typical  situation  to  be 
considered  in  applications  to  coating  problems  as  shown  in  part  A of  Figure  1-1 -(a)  and 
part  B of  Figure  1-1 -(b).  Now  the  viscosity  of  the  lubricant  may  change  considerably 
across  the  thin  film  as  a result  of  temperature  variations  that  arise  in  some  rotating 
systems  such  as  the  constraint  withdrawal  problem.  And,  the  pressure  also  changes  along 
the  thin  film  due  to  convergent  or  convergent-divergent  geometry  as  shown  in  part  C of 
Figure  1-1 -(b).  Using  Equation  (2.14)  is  satisfactory  in  majority  of  fluid  film  applications 
but  is  not  sufficient  for  usage  in  the  case  of  severe  operational  condition,  and  complicated 
boundary  conditions.  Considering  the  core  zones  which  may  be  attached  to  both  surfaces 
in  the  system,  this  approach  is  not  useful.  This  study  will  be  dealt  with  Chapter  7.  With 
these  expressions.  Equations  (2.14)  and  (2.15),  for  the  velocity  gradients  and  the  velocity 
components,  can  be  now  derived  as  expressions  for  the  volume  rate. 

The  volume  rates  per  unit  width  in  the  x and  y directions  are  defined  as 

q^  = £udz  and  q^  = ^vdz 


(2.16) 
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Substituting  Equations  (2.14)  and  (2.15)  into  Equation  (2.16)  gives 


r* 

9.  = Jo 


\ \ dp  ^ 

- — +Ci^  + C2 

dz  and  q = \l 

\ \ dp 

^ ^ z'  + CxZ  + Ci 

2 p dx 

2 p dy 

dz  (2.17) 


The  component  of  these  equations  represents  the  well  known  Poiselle  (or  pressure)  flow, 


and  the  Coutte  (or  velocity)  flow. 


2.3  Reynolds  Equation 


The  introduction  of  a fluid  film  between  moving  surfaces  provides  a lubricating 
effect  which  substantially  reduces  the  coefficient  of  friction  between  them.  As  previously 
discussed,  such  films  are  usually  sufficiently  thin  that  viscous  forces  are  large  in 
comparison  with  inertia  forces,  so  that  the  latter  may  be  neglected.  The  Navier-Stokes 
equation  in  fluid  dynamics  can  be  shown  to  reduce  to  the  Reynolds  equation  in 
lubrication  theory  when  inertia  effects  are  negligible  [22]. 

Using  the  preceding  Equations  (2.10)  and  (2.11)  as  a start  and  integrating  them 

with  respect  to  z with  the  boundary  conditions  that  w(/i)  = v(/i)  = 0 and 
u(0)  = U,v{0)  = V gives 

u = - — — z(z  -h)  + U — - — 

2ili  cx  n 


V = ——z[z  -h)  + V — - — 
2ju  0 h 


(2.18) 


Now  assuming  p ^ p(z)  and  integrating  the  continuity  equation  (compressible,  unsteady 
state  flow)  with  respect  to  z across  the  passage. 


dpu  rhdpv 


dp 


Ig—^dz  + \g—~dz  + pw{h)  + ~h  = 0 


dc 


dy 


dt 


(2.19) 


where  mh)  = — , Equation(2.19)  can  then  be  rewritten  as: 

dt 


12 


f/,4^  , (h^P^  , ^ph  ^ 

in  ~i:rdz  + L ~iT'dz  + — — = 0 


dx. 


dy 


dt 


(2.20) 


Now  applying  Leibnitz’s  rule  for  differentiating  an  integral  gives: 


= f[b{t)A^  - fWA^  + 


dx 


(2.21) 


The  result  of  this  operation  gives 


^ f/,  , ^ f*  j 

-ip^z  + -tpvd.  + — = 0 


(2.22) 


Since  p ^ p{z)  is  assumed,  it  can  be  removed  from  the  integral,  and  from  Equation 


(2.18) 


[h  j ^ .u  I ^ ,Vh 

igudz  = -——h  +— and  \oVdz  = ———h  + — 
12  ju  ac  2 12  n oy  2 


(2.23) 


Substituting  (2.23)  into  (2.22)  the  result  is  as  follows  [34], 


dx. 


ph^  d ( ph^\  4^ 

— + — — = 6 

\ H ) dx  dy\  ^ ) dy 


+ 72 


dph 

dt 


(2.24) 


The  physical  significance  of  equation  (2.24)  is  that  the  generation  of  hydrodynamic 
pressure  within  the  film  is  a function  of  three  terms  which  have  been  called  the  wedge, 
stretch,  squeeze  contributions  to  load  and  local  expansion.  The  wedge  term  appears  to  be 
the  most  important  of  the  three  terms  in  common  application  such  as  bearing  design. 
Most  lubrication  examples  assume  smooth  and  rigid  bounding  surfaces  for  the  lubricant, 
and  in  such  cases  the  stretch  term  has  no  significance.  The  squeeze  contribution  is  due  to 
normal  approach  of  the  surfaces,  whether  or  not  sliding  motion  occurs  between  them. 
This  effect  can  be  caused  by  slow  speed  squeezing  action,  or  by  relative  normal  vibration 


between  the  surfaces. 
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2.4  Energy  Equation 


During  the  lubrication  process  energy  is  dissipated  in  the  two  ways,  the  action  of 
viscous  friction  and  compression  (smaller  effect)  [35,  36].  Both  of  these  sources  release 
heat  which  changes  the  fluid  temperature.  Energy  dissipation  is  highly  dependent  on  the 
viscosity  of  the  fluid  which  changes  the  dissipation  as  well  as  the  flow  and  pressure.  The 
energy  equation  results  from  a balance  between  the  heat  generated  by  mechanical 
dissipation  and  the  heat  transported.  The  mechanisms  of  heat  transport  are  conduction 
and  convection.  The  effects  of  heat  conduction  are  described  by  means  of  Fourier’s  Law. 

A fundamental  equation  of  hydrodynamics  is  the  law  expressing  the  conservation 
of  energy  at  a point  in  the  lubricant  film.  In  the  unsteady-state  form,  with  no  external 
sources  of  film  heating  and  constant  pressure  across  the  film  thickness  it  may  be 
expressed  in  the  Cartesian  coordinates  as  [35,  37]: 


P 


PE  dE  PE  PE 
Pc^^  P^^^  Pz 


+ PP 


u 


Px 


A 


+ V 


7' 


+ w 


Pz 


-1 

\p) 


p_ 

Px, 


Pk)  Py 


f 


y A) 


+ 


(3z 


r 


\ 


.4; 

V Pz) 


X 


Pu  Pv  Pw 

Px.  py  Pz 


+ 0 


(2.25) 


The  first  term  is  the  change  in  the  internal  energy  and  internal  energy,  E,  which  can  be 
expressed  in  terms  of  known  quantities  of  lubricant.  Enthalpy  per  unit  mass,  V , is 


defined  as 


1 

1=E+V=E+- 

P 


and,  differentiating,  PE  - PI  - pW  - VPp 


Also 


pPV 


Define  the  coefficient  of  thermal  expansion  a as 
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a - 


V 


\jr) 


Therefore: 


dE  = CpT-aTVcp-p3^ 


also 


^ = CpT  - aT  4?/ p - p^lj  p) 


(2.26) 


Substituting  for  in  Equation  (2.25), 


pcp 


JT  JT  JT 

a: 


aT 


4?  4*  4^ 

M — + V— + W — 

3c  ^ a 


<3c 


^ a 
+ 


ar 


\ 


\ dx.)  4^/ 


+ 


a(  ar) 
a\  a) 


-A 


a aif  aw 

a 


+ 0 (2.27) 


where  the  last  two  terms  are  the  dissipation  function.  This  is  given  as  the  rate  at  which 
the  viscous  forces  do  work.  A is  volume  viscosity  and 


0 = p 


^ aiT 


+ 2 


f-T 


+ 2 


^ a-> 


2 

3 


\ 


a a aw 
ya  4 ^ 


+ 


a a 

— + — 

4 a 


\2 
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a aw 

— + — 

a 4 


\2 
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aw  a) 

— + — 

a a J 


The  physical  significance  of  equation  (2.27)  is  that  the  generation  of  temperature 
within  the  film  is  a function  of  four  parts  which  are  the  compression,  viscous,  convection, 
and  conduction.  The  first  term  in  the  left  hand  side  of  Equation  (2.27)  expresses  the 
convection  of  the  internal  energy  of  the  fluid.  The  convection  cooling  term  is  generally 
not  considered  to  be  significant  in  such  a lubrication  problem  with  very  narrow  passage. 
Cameron  justified  its  principal  cooling  term  with  an  order  of  magnitude  ^lnalysis  by 
comparing  it  in  application  to  a hydrodynamic  type  bearing  [34].  Accordingly,  the 
conduction  cooling  term  in  the  first  term  of  the  right  hand  side  is  more  important.  The 
relative  importance  of  the  heating  terms  may  provide  a useful  tool  to  avoid  some 


difficulty  in  solving  the  equation.  Results  estimated  with  the  viscous  heating  term  alone 
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as  the  heating  source  are  found  in  many  references  [35-37].  However,  in  this  study  all 
terms  are  considered. 


2.5  Summary  and  Conclusion 

This  chapter  has  dealt  with  the  basic  lubrication  equations  used  in  the  study  of 
thermohydrodynamic  problem.  The  Navier-Stokes  equations  have  been  derived  in 
general  form  for  lubricating  films.  The  Reynolds  equation  was  then  derived  from  the 
coupling  of  the  reduced  general  form  with  the  continuity  equation.  The  Reynolds 
equation  contains  viscosity  and  density  as  parameters.  These  properties  are  depend  on 
the  pressure,  film  thickness,  and  time.  The  energy  equation  governs  the  temperature 
distribution  within  the  lubricant.  The  coupling  of  the  Reynolds  equation  and  energy 
equation  will  be  dealt  with  in  Chapter  5,  where  the  full  theory  of  thermohydrodynamic 
lubrication  is  developed.  Velocity  expression  derived  from  the  reduced  general  form  is 
used  for  determining  the  shear  zone  using  Newtonian  behavior.  This  will  be  classified 
with  two  different  approaches  considered  to  determine  the  shear  zone  in  the 
thermohydrodynamic  problem. 
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Figure  2. 1 Schmatic  of  free  surface  flow 


CHAPTER  3 

FREE  SURFACE  FILM  FLOWS 
3.1  Introduction 


The  first  theoretical  treatment  of  the  withdrawal  processes  was  presented  by 
Jeffreys  [6].  Since  he  neglected  the  effect  of  surface  tension,  his  results  predict  the 
maximum  film  thickness  for  a given  set  of  conditions.  Surface  tension  forces  tend  to 
squeeze  out  liquid  from  the  film  thus  decreasing  film  thickness. 

A more  accurate  description  of  the  withdrawal  process  was  presented  by  Landau 
and  Levich  [3].  The  initial  work  of  Landau  and  Levich  regarded  the  problem  as  one- 
dimensional, neglecting  the  velocities  in  the  meniscus  region  normal  to  the  moving  sheet. 
Furthermore,  the  gravity  forces  were  also.  They  introduced  the  static  meniscus  condition, 
considering  equilibrium  between  gravitational  and  capillary  forces  in  the  thick  meniscus 
region.  Therefore,  their  work  considered  the  effect  of  surface  tension  and  made  the 
assumption  that  the  convective  momentum  terms  were  negligible.  The  curvature  match 
method  used  by  Landau-Levich  and  others  to  derive  withdrawal  theories  is  based  on 
matching  surface  curvatures  in  the  meniscus  region  [3].  Later,  White  and  Tallmadge 
included  gravitational  forces,  while  retaining  the  one-dimensional  assumption  [5,  8]. 
However,  in  order  to  apply  the  static  meniscus  condition  derived  by  Landau  and  Levich 
and  obtain  an  analytical  but  approximate  expression  for  the  relationship  between  the 
parameters  of  the  problem.  Soroka  and  Tallmadge  presented  the  solution  including  the 
nonlinear  inertia  effects  for  a Newtonian  fluid  [38].  They  assumed  that  the  non-linear 
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inertial  terms  in  momentum  equations  can  be  linearized  by  considering  the  treatment  of 
Oseen  [33].  It  is  found  that  the  inertia  effects  become  very  important  for  capillary 
number  in  the  case  of  Ca  > 1.0  [13-15]. 

The  previous  predictions  for  the  film  thickness  in  the  case  of  angular  withdrawal 
of  flat  plate  from  liquid  container  were  based  on  a one-dimensional  linear  approach 
introduced  by  Deryagin  and  Levi  and  developed  by  Tallmadge  [7,  8].  Deryagin  and  Levi 
obtained  a solution  by  assuming  constant  surface  tension,  a one-dimensional  flow  in  the 
thin  meniscus,  and  a negligible  value  of  the  thin  thickness  gradient  dh/dx«l  [9, 10]. 

There  have  been  many  studies  of  the  film  thickness  on  a rotating  cylinder. 
Tharmalingam  and  Wilkinson  applied  the  theory  of  flat  plate  withdrawal  to  that  of  the 
cylinder.  They  also  analyzed  the  thickness  of  the  liquid  adhering  to  cylinder  although 
their  predicts  are  too  large  [10]. 

When  a partially  submerged  cylinder  rotates  in  an  oil  bath  as  in  Figure  1.1 -(a)  a 
thin  film  of  liquid  is  withdrawn  from  the  container.  This  simple  system  provides  a model 
both  physically  and  mathematically.  It  is  considered  the  phenomenon  of  fluid  dynamics 
under  relatively  idealized  conditions.  It  is  desired  to  predict  the  flow  rate,  q,  which  is 
constant  around  the  periphery,  and  the  film  thickness,  h,  which  depends  on  the  angle  of 
measure,  dm,  as  functions  of  the  fluid  properties,  /u,  p,  and  cr  the  angle  of  contact  Oc,  the 

cylinder  radius,  R,  and  the  rotating  speed,  U.  The  withdrawal  problem  of  Newtonian 
liquids  onto  a rotating  cylinder  reduces  to  the  drag-out  problem  due  to  a plate  inclined  at 
an  appropriate  angle  Oq  to  the  horizontal  liquid  surface.  The  problem  for  the  inclined  flat 

plate  is  sketched  in  Figure  (3.2). 

The  study  described  in  this  proposal  deals  with  the  dynamics  of  liquid  films  on 
flat  plate  and  extends  the  case  to  the  cylindrical  geometry, 

1 . The  fluid  is  incompressible  and  isotropic 

2.  The  flow  is  laminar  and  stable 

3 . Gravity  is  the  only  body  force 
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4.  The  solid  support  is  completely  wetted  by  the  fluid 

5.  The  fluid  is  Newtonian 

With  contact  angles  other  than  zero,  and  using  the  assumptions  stated  above  the 
steady  state  momentum  Equations  (2.11),  (2.12),  and  continuity  equation  (2.13)  referred 
in  Chapter  2 becomes. 


u — + w — 


1 ^ 
p dx 


1 ^Xxz 

p dz 


- g sin 


(3.1) 


dz 


= pg  cos  P, 


(3.2) 


du  dv  ^ 
— + — = 0 

dx  dz 


(3.3) 


It  is  assumed  that  the  interfacial  shear  stresses  are  negligible  and  that  the  only  pressure 


existing  at  the  interface  besides  atmospheric  is  that  due  to  surface  tension: 


f -S  = Txz  = Q P = Pa 


atz  = h 


(3.4) 


The  free  surface  can  be  described  by  the  h(x),  and  the  radius  of  curvature  of  the  surface  is 


given  by  — = - t 


h!  dx‘ 


R 


l+{dh/dx)‘ 


!/2 


[3].  Since  the  curvature  is  assumed  to  be  small. 


Laplace's  capillary  statics  equation  can  be  used  to  approximate  the  pressure  [3,4]. 

(fh 

= (3.5) 

dx 


For  another  condition  at  the  free  surface  it  is  also  assumed  that  the  interfacial  shear 
stresses  exists  and  that  the  pressure  existing  at  the  interface  besides  atmospheric  is  that 
due  to  surface  tension  and  normal  stress: 

f -n  = rzz  = -r«  and + = atz  = h (3.6) 

where  r -s  and  r -n  are  the  shear  and  normal  components  of  stress,  respectively. 

At  the  wall  of  the  moving  plate  the  non-slip  condition  is  involved: 

u = U and  w = 0 atz  = 0 (3.7) 

The  Newtonian  relationship  between  shear  stress  and  rate  of  deformation: 
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Txz  = 


du 

dz 


(3.8) 


Combining  Equations  (3.4),  (3.5),  (3.6),  (3.8)  and  substituting  into  Equation  (3.1)  gives: 


du  du  <7  (fh  u d^u 

+ = — 7T  + — ^-gsin^, 

dx  dz  p dx  p dz 


(3.9) 


3.2  Landau-Levich  Solution 


In  the  theoretical  development,  the  withdrawal  process  is  divided  into  three 
regions  as  described  by  Gutfmger  and  Tallmadge,  Landau-Levich  [2,  3].  These  three 
regions  describe  different  phenomena  of  flow  dynamics  as  follows: 

Region  1 . This  is  the  region  above  the  liquid  level  where  the  film  thickness  is  a 
constant  (hoy  Surface  tension  effects  are  not  effective.  The  resulting  flow  equation  is: 


p d^u 
p dz!' 


-gsin6>c  = 0 


(3.10) 


Region  2.  This  is  the  "dynamic  meniscus"  region.  Here  viscous,  surface  tension, 
and  gravitational  forces  are  present: 


p d^u  cr  d^h 
p dz!  p dx 


gsin6>c  = 0 


(3.11) 


Region  3.  This  is  the  "static  meniscus"  region  near  the  liquid  surface.  Only 
surface  tension  forces  are  considered.  The  process  is  described  by  Laplace's  equation  of 
capillary  statics: 


cr 


cr 


(3.12) 
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An  approximate  solution  to  the  continuous  withdrawal  process  results  from  a 
combined  solution  of  the  three  regions.  Boundary  conditions  for  the  moving  film  are 
given  as 

= 0 atz  = h 

u-U  at  z = 0 (3.13) 

Equation  (3.10)  using  the  boundary  conditions  gives  the  velocity  profile  for  region  1 and 
the  flow  rate  can  be  obtained  by  integrating  Equation  (3.9). 

In  similar  manner  Equation  (3.11)  gives  the  flow  rate  for  region.  These  solutions 
to  regions  1 and  2 are  matched  by  using  the  steady  state  condition  of  equal  flow  rates  in 
both  regions. 


Uho- 


pg  sin  6>, 
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cr  d^h  pgsme, 
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yv  3 y 


(3.14) 


Equation  (3.14)  is  a third  order  non-linear  differential  equation  in  h with  the 


boundary  conditions, 


h = ho,  dh/dx  = 0,  h/dx^  = 0 


atx  = 0 


(3.15) 


This  equation  is  put  into  dimensionless  form  by  considering  the  treatment  of  White  with 
the  substitutions  [5], 


h 

$■  = — , To  = ho 
ho 


/ . \l/2 

pg  sm0c 
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pU 


yl 


,C  = ^,l  = -f(3Ca> 

O'  ho 


/3 


(3.16) 


With  the  substitutions.  Equation  (3.14)  becomes  in  dimensionless  form. 


dX^  3 ’ 


(3.17) 


with  the  boundary  conditions. 


dg 

dX 


= 0 


d^g 

dX^ 


= 0, 


at  X = 0 


(3.18) 


To  match  the  solution  for  the  static  meniscus,  region  3,  with  that  of  the  dynamic 
meniscus,  region  2,  the  condition  must  be  needed  that,  the  curvature  computed  from  the 
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dynamic  meniscus  Equation  (3 . 1 1 ) is  identical  to  that  computed  from  the  static  meniscus 
Equation  (3.12).  It  is  approached  from  the  dynamic  meniscus  side  or  from  that  of  the 
static  meniscus. 

The  curvature  match  condition  is  stated  in  dimensionless  form  as  [3,  5,  8], 


dA^  j 


^ d'^g^ 


dA, 


J r 


(3.19) 


Integration  of  the  Equation  (3.12)  for  the  region  3 with  the  boundary  condition  h = oo  at 
the  horizontal  liquid  surface  gives  [39], 


dh 

dx 


dh^ 
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(N 

1+ 

ydx  ^ 

-|l/2 
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(l-cos6>c)-l 


(3.20) 


From  the  condition  that  dh/dx  ->  0 as  h —>ho,  the  solution  of  Equation  (3.20)  is 


X — > 
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\l/2 


V 


/^(l-cos6>c) 


. Substituting  this  solution  into  Equation  (3.12),  the  curvature 


y 


is  obtained  as  follows, 
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dx' 


-> 
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which  becomes  in  dimensionless  form. 
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In  the  Equation  (3.17)  with  To^  <<  1 it  takes  the  form  as  follows, 

id^g  , 


(3.23) 


From  a numerical  integration  of  Equation  (3.23)  results  are  shown  in  Figure  (3.3)  as 
follows  [3,  5,  9], 


23 


j 
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Substitution  of  Equations  (3.22)  and  (3.24)  into  Equation  (3.19)  gives, 
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Introducing  the  dimensionless  variables,  the  Landau-Levich'  result  is, 
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(3.24) 


(3.25) 


(3.26) 


3.3  Extension  to  the  Landau-Levich  Model 


In  deriving  Equation  (3 . 1 7),  Landau  and  Levich  neglected  the  drainage  down  the 
plate  due  to  gravitational  force  [3].  When  gravity  is  taken  into  account,  the  governing 
equation  and  boundary  conditions  are  already  obtained  in  Equations  (3.17)  and  (3.18)  in 
dimensionless  form. 


To  match  boundary  conditions  of  Equation  (3.17)  when 


d^g 

dX^ 


1 2 


as  g 00, 


Equation  (3.17)  can  be  linearized  by  substituting  g = 1 + e and  neglecting  terms  of  the 
order  s [5].  Therefore, 


(3.27) 


After  making  a change  of  variables  to  X = 


3d^g 

■ yr 
dx 


= s 


, Equation  (3.27)  transforms  to 

(3.28) 


which  is  similar  to  the  Landau  and  Levich  equation.  In  same  manner  used  in  Equation 
(3.23), 
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(3.29) 
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where  S is  constant  and  independent  on  Ca  and  Tq- 

Equation  (3.29)  also  can  be  written  by  a change  of  variable  used  above  as, 


d^g 


dX 


/3 


1 


(3.30) 


Proceeding  using  the  same  method  as  Landau  and  Levich,  an  equation  that  related  the 
dimensionless  film  thickness  to  the  capillary  number  can  be  obtained.  Substituting 
Equations  (3.29)  and  (3.30)  into  Equation  (3.19), 
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The  modified  treatment  is  presented  to  give  a more  realistic  boundary  condition  [11,  12]. 
The  tangential  stress  at  the  free  surface  becomes. 


^ {txx  - Tzz)sin  20  + cos  20  = 0 


(3.32) 


and  the  normal  stress  at  the  free  surface  becomes. 


P + Txx  sin^  0 + Tzz  cos''  0 + Txz  sin  20  = p^~ 


a 


(3.33) 


The  boundary  conditions  apply  in  the  case  of  0=  (f  due  to  the  free  wave,  thus  Equations 
(3.32)  and  (3.33)  become  as  follows. 


^ o 

Txz  = -M—  = 0 

oz 
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To  give  the  velocity  profile  from  equation  (3.1 1)  as. 


(3.35) 
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Differentiating  Equation  (3.35)  and  substituting  Equation  (3.35)  into  Equation  (3.36),  the 
result  is  given  by 
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dx^  dx 
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Equation  (3.38)  can  be  linearized  by  substituting  g=l  + s where  s is  small  and 
neglecting  terms  of  the  order  [5].  Introducing  the  change  of  variable  again,  Equation 
(3.38)  becomes. 
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Upon  numerical  solution  of  Equation  (3.39)  the  curvature  at  j ->  0 has  a finite  value. 


f \ 

d^g 


dx 


= S{Ca,T,) 


Jx^o 


(3.40) 


The  numerical  integration  of  Equation  (3.40)  can  be  carried  out  using  Euler 
predictor-corrector  method  with  the  starting  value  procedure  introduced  by  Landau  and 
Levich  [5].  Values  oiS(Ca,  Tq)  are  plotted  in  Figure  (3.2)  with  Tq  and  Ca.  Here,  at  the 

high  Capillary  numbers,  the  results  show  that  film  thickness  are  relatively  larger  as  the 
values  of  the  curvature  coefficient  S rapidly  increase. 

Substituting  Equation  (3.40)  into  Equation  (3.22),  therefore. 
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The  above  Equations  (3.42)  have  been  tested  over  a range  speed  and  viscosity 
indicated  by  a range  of  Capillary  Number  from  0.001  to  10  as  given  in  Table  2.  There  is 
no  reason  to  doubt  its  applicability  at  value  of  Ca  lower  than  0.01.  Calculated  data 
obtained  from  three  theories  are  considered  to  be  good  when  0.01  < Ca  < 0.05  as  shown 
in  Figure  (3.5).  The  small  differences  may  not  be  significant  since  they  are  not  affected 
by  any  surface  tension  and  gravity  force  in  this  region.  For  high  Ca,  from  0.05  to  0.5  the 
results  of  Equation  (3.42)  with  realistic  boundary  conditions  seems  to  have  more  accurate 
values  compared  with  the  other  two  results  from  Equations  (3.26)  and  (3.31).  However, 
for  higher  Capillary  Number,  Ca  > 0.5,  the  results  of  Equation  (3.36)  seem  to  have  stable 
values  compared  with  other  two  results  from  Equations  (3.31)  and  (3.42).  Significant 
differences  in  this  region  are  due  to  the  assumption  of  negligible  gravitational  forces  in 
Equation  (3.26)  with  respect  to  viscous  forces,  and  inertia  forces  in  all  theories.  Thus, 
when  Ca  > 0.5  it  is  necessary  that  the  theoretical  development  include  inertia  forces  to 
obtain  results  much  better  than  that  given  by  previous  theories.  Some  assumptions  may 
be  limited  above  much  high  Capillary  Number  by  the  occurrence  of  waves  or  turbulence 


even  though  inertia  forces  are  considered. 


Table  2:  Constant  Values  Dependent  on  Ca  and  To 


Ca 

H 

o 

U) 

“S(CaJop^ 

Ca 

0.001 

0.1954 

0.6446 

"■■'0.6 

0.003 

0.2330 

0.6473 

1.0 

0.006 

0.2602 

0.6503 

zo 

0.01 

0.2823 

0.6536 

3.6 

0.03 

0.3369 

0.6660 

4.0 

0.06 

0.3775 

0.6797 

6.0 

0.1 

0.4113  ' 

0.6945 

8.0 

0.2 

0.4635  ' 

0.7240 

10.0 

0.3 

0.4979 

0.7481 

Note  : (3)  Equation  (3.42) 


1 To  (3) 

S(Ca,To) 

0.5638 

0.8054 

0.6177 

0.8651 

0.6956 

0.9791 

0.7417 

1 .0694 

0.7736 

1.1469 

0.8163 

1.2798 

0.8443 

1.3944 

0.8644 

1.4971 
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3 .4  Application  to  the  Cylindrical  Coordinates 


Equations  (3.26),  (3.31),  and  (3.42)  have  the  non-dimensional  uniform  thickness. 
To,  to  be  calculated  as  a function  of  the  Capillary  number,  Ca,  and  the  angle  of  contact 

Oq.  By  using  the  final  uniform  film  thickness,  ho,  the  liquid  flow  rate,  q,  and  the  steady 

state  condition  of  equal  flow  rates  can  be  calculated  as  follows. 


q = Uho- 


pgho  sin6>c 


(3.43) 


Knowing  the  predicted  flow  rate  based  on  the  uniform  film  thickness  it  is  possible  to 
apply  it  to  the  rotating  cylinder  for  the  corresponding  contact  angle.  The  analysis  which 
is  applicable  to  the  flat  plate  of  Equations  (3.26),  (3.31),  and  (3.42)  can  be  modified  to 
determine  the  positions  and  the  variation  of  film  thickness  around  the  cylinder  as  follows 

[9], 


q = Uh 


pgh^sin0„, 

3/j 


(3.44) 


where  h is  the  film  thickness  at  any  point,  0^  and  R>>h. 


At  the  top  of  the  rotating  cylinder,  Om  = 0,  Equation  (3 .44)  gives. 


q = Uh  (3.45) 

where  h is  the  film  thickness  on  the  top  of  the  rotating  cylinder. 

Using  Equation  (3.43)  throughout  Equation  (3.45),  it  gets 


h^ho- 


pgh  ^s\n0, 
3//f7 


pgh^  sin  0„, 
3juU 


(3.46) 


Multiplying  Equation  (3.46)  by  (pgsin0c/ju)''^U,  an  alternative  form  is  as  follows. 


T=To- 


(3.47) 
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where 


T = h 


pg^vciOc 


pU 


J 


Similarly,  multiplying  Equation  (3.47)  by  (pgsinOc/p)^'^U,  the  form  is, 
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T =T — T 
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sin  e„, 


\ 


sin^ 


(3.48) 


c J 


where  T-h(pgsinOJ is  the  film  thickness  at  the  position  on  the  rotating  cylinder. 
The  dimensionless  film  thickness,  T,  can  be  calculated  by  using  the  Newton-Raphson 


method. 


3.5  Mathematical  Solution  and  Discussion 


There  is  experimental  data  with  which  the  final  results  of  the  Equation  (3.48) 
given  by  Equations  (3.26),  (3.31),  and  (3.42)  can  be  compared.  That  is  the  direct 
measurement  of  steady  state  thickness  of  a film  withdrawn  on  a rotating  cylinder  as 
carried  out  at  the  University  of  Wisconsin-Madison.  The  oil  used,  its  properties,  and  the 
measured  thickness  are  given  in  Table  3.  This  is  accomplished  by  scraping  fluid  off  a 
continuously  rotating  cylinder  as  shown  in  Figure  (3.1).  This  experimental  device 
consists  of  scraper  and  rotating  steel  disc.  Oil  attached  to  the  surface  of  rotating  disk  is 
collected  by  scraper  located  at  the  side  of  the  disk.  Theoretical  values  of  dimensionless 
thickness,  T,  for  various  Capillary  numbers,  Ca,  are  changed  to  the  dimensional  form  for 
comparison  with  the  experimental  results. 


Table  3:  Experimental  Data  and  Lubricant  Properties 


Motor  Oil 
(23.90C) 

rpm 

h (cm) 
(6m=  -900) 

g 

(m/s^) 

(N-sW) 

P 

(kg/m^) 

a 

(N/m) 

SAE  low 

40 

0.08382 

9.807 

0.08 

878.24 

3.6x10-2 

90 

0.13462 

180 

0.17780 

Source:  (personal  communication,  A.  A.  Seireg,  1996) 
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The  calculated  results  by  three  different  theories  are  presented  graphically  in 
Figures  (3.6)  to  (3.8)  by  plotting  the  thickness,  h,  against  rotating  speed,  rpm.  As 
expected  the  results  show  that  the  film  thickness  becomes  thinner  as  the  position  angle, 
6ni  decreases.  It  is  considered  that  the  results  are  more  reliable  for  values  of  Gc  of  90® 
due  to  the  characteristics  of  the  assumed  model.  Hence,  the  numerical  solution  is  not 
stable  beyond  ^^=30°  for  equations  (3.26)  and  (3.42).  As  rotating  speed,  N (rpm), 
increases,  the  film  thickness  also  increases  in  the  low  speed  range.  The  dimensionless 
film  thickness  depends  on  one  dimensionless  physical  parameter  Ca.  When  the  surface 
tension  of  the  lubricant  is  given,  the  value  of  Ca  is  related  to  the  value  of  the  product  juU. 
As  the  velocity  increases  or  decreases  Ca  corresponds  to  a variation  of  meniscus  shape  at 
low  speed.  The  meniscus  shape  for  various  speeds  changes  the  film  thickness  attached 
on  the  rotating  cylinder  surface. 

Figure  (3.9)  shows  a comparison  of  these  theories  with  the  experimental  data 
given  in  Table  3 over  the  range  of  the  experimental  data  and  the  theoretical  results  listed 
in  Table  2.  It  is  seen  that  the  results  of  Equation  (3.31)  agree  very  closely  experimental 
data  up  to  a value  of  about  100  rpm  (Ca=l  .7733  in  dimensionless  value).  The  differences 
between  the  theories  and  experimental  data  at  high  speed  (or  Ca)  is  to  be  expected  due  to 
the  assumption  of  one-dimensional  flow  and  negligible  inertia  terms  in  the  meniscus 
region.  For  higher  values  of  Ca  such  assumption  is  no  longer  valid  and  capillary  statics 
are  no  longer  adequately  describe  the  shape  of  the  film  near  the  liquid  container. 
Accordingly  in  this  study  the  theoretical  prediction  of  Tallmadge  and  White,  Equation 
(3.36),  will  be  mainly  considered  to  analysis  the  effect  of  surface  tension  in  the 
phenomenon  of  the  lubrication  problem  between  two  surfaces  system  which  will  be 
discussed  later. 

Figures  (3.10)  and  (3.11)  show  the  film  thickness,  h,  at  different  circumferential 
positions  on  the  rotating  cylinder,  0^,  for  the  SAE  lOW  motor  oil.  It  can  be  seen  in 

Figure  (3.10)  that  the  larger  rotating  speed  of  the  cylinder,  U,  the  larger  the  film 
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thickness,  h,  is.  Figure  (3.10)  also  shows  that  the  film  thickness,  h,  decreases  with 
decreasing  6m  for  contact  angle,  6c  = 1 20®.  This  can  be  attributed  to  the  increase  in  the 
velocity  of  the  film  in  the  circumferential  direction  due  to  gravity  forces.  The  contact 
angles,  6c  = 60®,  90®,  and  120®,  which  correspond  to  the  situations  of  the  cylinder 
submerged  to  three  different  levels  are  considered.  Figure  (3.1 1)  shows  the  variations  of 
film  thiekness,  h,  against  6m  for  these  contact  angles.  The  film  thickness,  h,  ean  be  seen 

to  approach  a constant  value  with  respeet  to  6m  near  the  top  of  the  cylinder. 

3.6  Summary  and  Conelusion 

A general  solution  of  the  continuous  withdrawal  process  from  the  liquid  container 
can  be  made  using  marching  technique  for  integrating  the  Navier-Stokes  equations  (one 
continuity  and  two  momentum  equations).  Some  approximations  and  simplification  of 
the  two  momentum  equations  are  introduced  based  on  physical  considerations.  A 
withdrawal  analysis  on  a rotating  eylinder  is  introduced  as  an  extension  of  the  flat  plate 
theory  of  Landau-Levich  without  the  effect  of  inertia  forees.  The  results  will  be  useful  in 
determining  lubrieating  zone  in  the  hydrodynamic  analysis. 

The  withdrawal  theory  is  not  entirely  satisfactory  for  predicting  the  film  thiekness 
at  high  speed  (N>100  rpm).  Numerical  analysis  by  the  Landau-Levich  method  is 
restricted  to  Ca<10  for  its  application.  It  is  apparent  that  this  is  due  to  the  reduced 
Navier-Stokes  equations  under  the  assumption  of  inertia  free  and  one-dimensional  flow. 
A more  detailed  study  of  the  two-dimensional  flow  with  the  effect  of  inertia  forces  and 
interfacial  boundary  conditions  is  necessary  to  improve  the  solution  for  much  higher 
speeds. 
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Figure  3.1  Diagram  of  a partially  submerged  rotating  cylinder  (disk) 

in  an  oil  batrh 
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(viscous,  surface  tension) 
Region  1 Region  2 Regions 


Figure  3.2  The  analogy  for  flat  plate  with  contact  angle  6c 


derivatives  of  ^ 
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Figure  3.3  Computed  solution  to  Equation  (3.23) 
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Figure  3.4  The  curvature  coefficient 


S is  Equation  (3.40)  as  function  of  and  Ca 
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Figure  3.5  Comparison  of  theoretical  predictions 

(constant  film  thickness  in  the  case  of 
a inclined  flat  plate) 
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Figure  3.6  Film  thickness  for  Landau-Levich  solution 

to  Equations  (3.26)  and  (3.48) 

0^=  131.81” 


h (inch) 
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Figure  3.7  Film  thickness  for  Equations  (3.31)  and  (3.48) 
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Figure  3.8  Film  thickness  for  Equations  (3.42)  and  (3.48) 
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Figure  3.9  Film  thickness  at  9^=  -90°,  0 =131.81° 

ill  ^ 

Comparison  between  analytical  results  of  the 
3 different  methods,  and  experimental  data 
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Figure  3.10  Film  thickness  for  Equation  (3.31) 

(Effects  of  rotating  speed  of  cylinder) 


h (inch) 


Figure  3.1 1 Film  thickness  for  Equation  (3.31) 

(Effects  of  contact  angle) 


CHAPTER  4 

FREE  SURFACE  FILM  FLOWS  IN  THE  PRESENCE  OF 

A STATIONARY  CONSTRAINT 

4.1  Introduction 

There  are  several  basic  problems  in  the  study  of  slow,  and  two-dimensional  flow 
of  a Newtonian  viscous  fluid  in  which  a single  stream  separates  when  constrained.  An 
example  is  shown  in  Figure  (4.1).  Here,  a rotating  disc  (or  cylinder)  and  a stationary  flat 
plate  are  partially  immersed  in  a large  container  of  viscous  fluid.  As  the  cylinder  turns, 
fluid  is  drawn  on  the  outside  of  rotating  cylinder  and  enters  the  narrow  gap  between  the 
two  bodies.  After  flowing  through  the  section  of  the  closed  area  between  the  cylinder 
and  flat  plate  the  fluid  separates  under  a surface  tension  membrane  into  two  films. 

The  system  in  the  above  case  has  a convergent-divergent  geometry  in  the 
direction  of  motion  and  the  film  adhering  to  the  moving  surface  will  be  dragged  into  the 
narrowing  clearance  space,  thus  building-up  a pressure  which  is  sufficient  to  carry  load. 
This  is  the  principle  of  hydrodynamic  lubrication,  a mechanism  which  is  essential  to  the 
efficient  functioning  of  many  systems  in  industry. 

The  problem  is  to  how  analyze  the  phenomenon  of  the  flow  under  the  surface 
tension  membrane  at  the  inlet  and  outlet  regions.  The  shape  of  the  membrane  is  not 
known,  but  must  be  found  as  part  of  the  solution.  The  mathematical  formulation  must 
consider  the  cylinder  rotation  speed,  the  viscosity  and  surface  tension  of  the  fluid, 
cylinder  radius,  and  clearance  between  two  surfaces. 
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4.2  Literature  Review 

The  withdrawal  process  by  two  surfaces  of  viscous  fluid  has  received  much 
attention  over  a period  of  half  a century.  The  reason  for  this  is  partly  due  to  the  difficulty 
of  defining  suitable  criteria  for  finding  the  position  of  the  line  of  separation  or  rupture  in 
an  adhered  layer  between  diverging  surfaces.  Current  interest  in  this  type  of  separation 
or  rupture  appears  to  have  started  with  an  investigation  concerning  the  hydrodynamic 
behavior  of  lightly  loaded  surface  in  a fluid-filled  capillary  container  [40-44].  A similar 
two-dimensional  separating  flow  at  the  exit  from  the  nip  of  an  eccentric  journal  bearing 
was  experimentally  investigated  by  Taylor  [43].  Approximate  theoretical  treatment  of 
separation  flows  were  undertaken  by  Pearson  [18],  and  Pitts  and  Greiller  [40].  In  both 
investigations,  gravitational  effects  were  not  included  as  their  influence  was  found  to  be 
small.  Some  analysis  for  these  flows  such  as  the  forward  rolling  process  where  the  liquid 
flowing  into  a narrow  gap  between  two  rotating  cylinders  have  dealt  with  the  problem 
from  three  perspective  (i)  classical  lubrication  theory  [21,  34,  40,  41],  (ii)  asymtotic 
expansion  [45],  and  (iii)  solution  of  the  Navier-Stokes  system  in  relevant  flow  domain  by 
Galerkin/fmite  element  method  [42].  For  fluid  film  systems  in  a convergent-divergent 
geometry  some  realistic  cavitation  boundary  condition  at  the  line  of  seperation  or  rupture 
have  been  suggested  by  Swift  and  Stieber,  Floberg,  and  Coyne  and  Elrod  [19,  46,  47]. 
They  are  applied  to  systems  such  as  the  rolling  system  in  viscous  flow  by  Taylor, 
Middleman,  and  Savage  [21,  41,  44].  In  this  chapter  the  withdrawal  process  of  the 
viscous  flow  in  the  cylinder-plane  geometry  is  modeled  using  lubrication  theory  with 
known  boundary  conditions.  The  prineipal  objective  is  to  extend  the  work  of  Middleman 
and  Savage  to  the  general  lubrication  situation  which  involves  flow  rupture  which  occurs 
at  the  exit  in  plane-cylinder  geometry  as  influenced  by  the  rotating  speed  and  surface 
tension  [21, 44]. 
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4.3  Mathematical  Formulation 


The  flow  pattern  in  the  region  of  the  meniscus  is  complicated  and  exact 
theoretical  treatment  of  the  problem  is  extremely  difficult.  The  motion  of  the  liquid  in 
the  region  where  the  gap  is  narrow  is  predominantly  determined  by  the  drag  due  to  the 
rotating  cylinder  and  the  variation  in  the  width  of  the  gap.  It  can  be  shown  that  inertia 
terms  in  the  hydrodynamic  equations  can  be  neglected  without  serious  error  except  at 
very  high  speed  as  illustrated  in  results  of  previous  chapter.  Terms  representing  the 
effect  of  gravity  are  included  in  the  original  treatment  but  can  be  also  omitted  for 
simplicity,  since  their  influence  is  small  except  for  the  lowest  speeds  of  rotation  and  for 
the  larger  gap  widths. 

To  calculate  the  flow  in  the  narrow  part  of  the  channel,  it  will  be  convenient  to 
use  the  coordinate  system  in  Figure  (4.1).  Introducing  the  dimensionless  variables  to 
lubrication  with  the  gravity  term  for  Newtonian  liquid  [16-18,  21], 
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y ^ h{x) 

T]  = , 0 = — — 
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Equations  (2.11)  and  (2.14)  of  the  two-dimensional  case  can  be  transformed  to  the 
following  forms. 


^ 


-r 


(4.2) 


d(D 

Equation  (4.2)  may  then  be  integrated  twice  with  respect  to  rj  to  give  u in  the  form 
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where  H(^)~  h{x)  / Ho  is  the  dimensionless  separation  between  the  two  surface.  In 
Equation  (4.4),  the  boundary  conditions  are  used  as  follows; 

</>  = 0 at  T]  = 0 


(4.5) 


It  is  necessary  to  introduce  an  additional  equation  for  the  conservation  of  mass.  This  is 


given  in  the  form 


q = U Hl  = f^^^^°udy 


Also,  a dimensionless  film  thickness  can  be  written  as 


0 = 


Hi 


Ho  UHo 


where  Hl  is  a constant  film  thickness  after  rupture  at  outlet  zone  and  then 


(4.6) 
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Hl  is  the  film  thickness  which  becomes  imiform  on  a region  of  the  cylinder  surface  far 

removed  from  the  meniscus  region  where  q is  the  flow-rate  entrained  by  the  cylinder. 
After  performing  the  integration  implied  by  Equation  (4.7)  and  rearranging  Equation 
(4.8)  it  can  be  solved  for  dP/d^  to  give: 


It  follows  then,  that 


dt] 
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Choosing  an  assumed  position  that  has  the  external  atmospheric  pressure  and  no 
gravity  force  at  the  outlet  region,  it  follows  then  that 

= o (4.12) 

and  then  Equation  (4.1 1)  is  reduced  to. 


P(^)  + r^  = i2jf  ? 


H 


(^) 


(4.13) 


where  h[^)  = 1 + — i^\  This  expresses  'no  slip'  on  the  cylinder  surface,  defined  by  the 


function  is  the  circular  arc 


(4.14) 


This  expression  introduces  algebraic  complications  and  a common  approximation  for 

, taking  note  of  the  fact  that  Ho/R«l  in  any  lubrication  problem  is  obtained  by 
expanding  the  bracketed  term. 

Equation  (4.13)  may  then  be  integrated,  and  the  result  is 
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0 is  still  an  undetermined  parameter.  To  obtain  a relationship  between  and  a 
boundary  condition  must  be  imposed  on  the  pressure.  From  Figure  (4.1)  it  can  be 
expressed  as  P{-^^  = 0 by  considering  a partially  flooded  inlet  boundary  condition  and 


calculated  from: 


0 = 12 
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Symbolically, 


0 = /(©,^,.#„r)  (4.17) 

Equation  (4.17)  expresses  the  relation  between  the  entry  position  and  the  rupture  position 
by  using  Swift-Stieber  boundary  condition  and  partially  flooded  condition  explained  as 
one  of  the  boundary  models  considered  in  the  next  section  (4.5). 


4.4  Coyne  and  Elrod  Boundary  Condition 

The  Swift-Stieber  or  Reynolds  boundary  condition  for  cavitation  applies  a zero 
pressure  gradient  between  the  full  film  and  the  cavitating  regions  [47].  It  is  known  to 
result  in  reasonably  precise  predictions  of  bearing  performance  for  moderately  or  heavily 
loaded  lubricated  contacts.  This  does  not  necessarily  imply  that  the  physical  assumptions 
of  the  condition  are  correct.  For  lightly  loaded  situations,  however,  the  boundary 
condition  proves  inadequate. 

The  occurrence  of  a reserve  flow  adjacent  to  a stationary  surface  causes 
separation  of  the  flow.  Coyne  and  Elrod  have  performed  a sophisticated  analysis  to 
determine  the  location  of  such  an  interface  and  interpreted  their  results  in  terms  of 
boundary  conditions  on  the  full  film.  Their  analysis  makes  the  following  assumptions 
[46], 

(1)  The  cavity  region  is  two-dimensional  and  corrugations  do  not  exist  at  the 
leading  edge.  All  the  lubricant  is  therefore  carried  away  by  the  surface  that  is  moving 
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and  separation  of  the  flow  occurs  from  the  opposite  stationary  surface.  The  analysis  is 
not  directly  applicable  to  the  situation  in  which  both  surfaces  are  in  motion. 

(2)  A balance  of  forces,  including  surface  tension,  exists  across  the  interface. 

(3)  The  interface  is  a free  surface. 

The  boimdary  condition  for  Reynolds  equation  considering  pressure  continuity  is 
expressed  as 


p = -—  + ^p  (4.18) 

r 

where  r is  the  radius  of  curvature  of  the  free  film  at  the  separation  point.  The  first  term  of 
Equation  (4. 1 8)  represents  the  pressure-drop  across  the  interface  due  to  surface  tension. 

Introducing  the  dimensionless  variables  into  Equation  (4.18)  the  boundary 
condition  represents  a force  or  pressure  balance  at  the  separation  point,  Xa , and  equates 
the  pressure,  p,  to  the  pressure  associated  with  surface  tension,  cr , 

(j 

p = — atx  = Xa  (4.19) 

r 


otherwise 


(4.20) 


where  the  additional  parameter  is  a modified  capillary  number  as  follows  [16-18], 
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(4.21) 


The  curvature  of  the  surface,  r,  is  unknown  a priori.  Since  the  surface  tension 
contribution  is  a small  one,  a simple  approximation  to  r proves  to  be  useful.  The  free 


surface  is  taken  to  be  semicircular,  then  a good  approximation  is 


r = h{xa)  - Hl 


(4.22) 


or 


r h{xg)  Hl 
Ho  ~ Ho  Ho 


(4.23) 


Substituting  Equation  (4.23)  into  Equation  (4.22),  we  obtain: 
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1 


at  ^ (4.24) 


4.5  Boundary  Conditions 
4.5.1  Partially  flooded  inlet  boundary 

The  fully  flooded  or  downward  inlet  boundary  for  large-sized  cylinder  (p  = 0 at 
Xt,=-oo)  in  Figure  (4.1)  cannot  occur  in  a real  system.  It  is  just  a convenient 
assumption  for  simplifying  the  analytical  solution  of  the  Reynolds  equation.  Most  of 
these  cases  are  presented  as  fully  flooded  for  known  positions  of  the  inlet  boundary.  In 
the  hydrodynamic  problems  this  assumption  has  widely  been  used  [18,  21, 42]. 

If  the  inlet  region  is  adequately  filled  with  the  lubricant  on  the  circumferential 
surface  entrained  by  the  rotating  cylinder  from  the  liquid  container  the  inlet  pressure 
buildup  will  act  on  the  free  surface  of  lubricant  beyond  a certain  point  such  that  the  inlet 
pressure  is  insensitive  to  further  increase  in  the  lubricant  supply  entrained  by  the  rotating 
cylinder. 

For  practical  cases,  the  extent  of  the  pressure  zone  is  small  compared  with  the 
dimensions  of  the  cylinder's  radius.  It  is  also  expected  that  the  pressure  does  not 
contribute  much  to  load  capacity  until  the  closest  zone  between  cylinder  and  the  flat 
plate.  For  a given  load,  speed  and  viscosity,  only  part  of  the  lubricant  injected  into  the 
nip  and  carried  on  the  cylinder's  entraining  surfaces  may  pass  through  to  outlet.  Such  a 
state  arises  when  there  is  a surplus  of  lubricant  in  entry  zone.  In  order  to  predict  the 
position  where  the  pressure  starts  to  develop  the  surface  tension  effects  are  considered  to 
be  negligible.  Accordingly  the  contact  point  where  the  free  surface  contacts  the 
stationary  flat  plate  becomes  the  starting  point  of  pressure  development. 

Pin  Palm 


at  X = - Xb 


(4.25) 


50 


4.5.2  Swift-Stieber  outlet  boundary  condition 

The  Swift-Stieber  (Reynolds)  outlet  boundary  condition  has  already  been  used  to 
solve  the  Reynolds  equation.  It  is  appropriate  here  to  discuss  it  in  more  detail  as  it  is  the 
most  frequently  employed  method  when  solving  this  problem.  In  the  experimental 
observations  that  prompted  its  formulation  for  relatively  high  speeds  and  thin  films,  there 
was  no  noticeable  negative  pressure  [16,  17].  The  outlet  was  marked  by  lines  of  shaped 
air  cavities  separated  by  liquid  streams.  The  cavities  are  considered  to  occupy  most  of 
the  clearance  space  depth  except  for  liquid  layers  on  the  moving.  Each  of  the  moving 
cylinder's  surface  leaves  with  a residue  of  evenly  spaced  circumferential  liquid  streamer 
of  lubricant.  The  appropriate  condition  at  the  film  rapture  can  be  formulated  by  flow 
continuity  considerations.  In  the  region  of  the  air  cavities,  the  pressure  is  taken  to  be  zero 
with  only  Couette  flow  (flow  caused  by  surface  velocities  only).  As  there  is  flow 
continuity,  a suitable  boundary  condition  for  the  general  case  is, 

/^  = — = — = 0 at  X = Xa  (4.26) 

ox  oy 

A boundary  condition  involving  surface  tension  would  have  to  be  introduced  for  more 
detailed  dynamic  analysis.  First  of  all,  it  is  not  necessary  to  formulate  a complex  model 
in  order  to  yield  results  for  the  general  features  of  this  system.  For  the  undetermined 
parameter,  0,  it  can  be  evaluated  from  the  additional  boundary  condition  on  pressure, 
imposed  at  ^ 


Accordingly, 


(4.27) 


0 = 12 


(4.28) 


Arranging  Equation  (4.28)  for  & , 


(4.29) 


Substituting  Equation  (4.29)  into  Equation  (4.13),  a relationship  is  obtained  for 

and  . In  the  absence  of  gravity,  the  values  of  and  can  be  calculated  by  using  a 
simple  numerical  technique. 

4.5.3  Separation  outlet  boundary  condition 

The  separation  type  outlet  boundary  condition  is  appropriate  at  low  speeds  and 
thick  films.  It  was  formulated  following  experimental  observations  of  a distinct  sub- 
cavity pressure  loop  prior  to  the  film  end  that  occurred  under  certain  conditions  [16,  17]. 
Thus,  it  was  considered  that  flow  separation  might  be  contributing  to  the  film  rupture. 
Separation  first  occurs  at  the  location  xa  of  the  outlet  film,  where  its  velocity  profile  has  a 

value 

u = —=0  (4.30) 

This  implies  that  the  position  is  both  a stagnation  point  and  a point  of  separation.  Such  a 
position  may  be  located  either  on  the  interface  or  at  the  onset  of  the  reverse  flow  region. 
The  pressure  gradient  is  positive  just  upstream  of  xa  with  the  pressure  having  a slight 

negative  value.  This  is  analogous  to  the  zero  reverse  flow  inlet  condition,  where  the 
reverse  flow  upstream  is  eliminated.  Here,  reverse  flow  downstream  in  the  expanding 
exit  causes  suspended  gas  in  the  lubricant  to  form  a cavity  pattern  of  closely  spaced  air 
fingers  more  rounded  at  their  ends  than  the  Swift-Stieber  type  cavities.  In  this  case,  the 
majority  of  the  lubricant  flow  is  believed  to  be  carried  in  layers  adjacent  to  the  entraining 
surfaces. 

For  simplicity,  Tl  = dP I -T , Equation  (4.2)  can  be  transformed  to  the 


following  form; 
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If  gravity  is  accounted  for  a trivial  modification  to  this  condition  subsequent  work  will  be 
done  in  same  manner  used  in  section  4.3.  The  velocity  is  given  as  follows. 


^(;;)  = H{^)t]]+  -J— 


(4.32) 


To  find  the  pressure,  it  is  necessary  to  satisfy  the  conservation  of  mass  of  Equation  (4.6). 
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(4.33) 


Changing  to  an  integration  form  the  pressure  gradient  is  found  to  be 


dn. 


\m)~® 


(4.34) 


The  boundary  conditions  on  pressure  asserts  that  the  pressure  including  gravity  effects  at 
inlet  region,  n , is  ambient  as  x ^ -Xt . A boundary  condition  at  the  outlet  region  is 


]_ 

n(^)  = 12j|^^ 


- 0 
h\^) 


(4.35) 


The  boundary  condition  represents  a force  or  pressure  balance  at  the  point  and 
equates  the  pressure  n to  the  associated  surface  tension  cr  as  shown  in  Equation  (4.20). 

When  the  boundary  condition  in  the  Equation  (4.20)  is  considered,  the  pressure 
balance  becomes  as  follows. 


1 


- 0 


(4.36) 


The  final  boundary  condition  expressing  stagnation  point  of  (j)  about  the  t]=0  and 
^ , and  expressing  a wetting  angle  of  PO"  is 

^ = 0 at  T]=0  and  4=  (4.37) 

orj 

When  boundary  condition  Equation  (4.37)  is  imposed,  the  results  of  Equation  (4.34)  is 
found  to  give 
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® = (4.38) 

To  obtain  a relationship  between  and  Coyne  and  Elrod  boundary  condition, 
Equation  (4.24),  must  be  imposed  on  the  pressure.  From  Figure  (4.1)  it  can  be  expressed 
as  P(-4b)  ~ 0 for  inlet  zone. 

4.6  Mathematical  Solution  and  Discussion 


A surplus  of  lubricant  in  the  entry  zone  makes  the  previous  model  difficult  for 
predicting  the  exact  position  of  the  inlet  boundary  due  to  the  fluid  mechanism  in  the 
convergent  geometry.  However,  if  there  exists  some  means  of  establishing  the  outlet 

boundary  distance,  and  inlet  boundary  distance,  from  the  closest  nip  in  the 

lubricating  zone  then  Equation  (4.16)  would  give  a solution  for  ©(f).  For  fluid  film 
plane  flat-cylinder  which  have  a convergent  clearance  space  geometry,  two  realistic 
cavitation  boundary  conditions  have  been  suggested.  In  the  Swift-Stieber  condition,  the 
cavity  pressure  usually  taken  to  be  atmospheric  and  surface  tension  effects  are  usually  of 
secondary  importance.  The  use  of  this  boundary  condition  in  analysis  predicts  pressure 
profiles  with  relative  ease.  Although  separation  boundary  condition  has  often  been 
suggested  for  lightly  loaded  fluid  film  bearings,  not  a great  deal  of  numerical  design 
information  appears  in  the  bearing  situation  to  be  discussed  in  the  next  chapter. 

Figure  (4.2)  shows  the  solution  for  ©(f)  of  Equation  (4.16)  obtained  by  using  the 
Swift-Stieber  cavitation  boundary  condition.  As  expected  the  effect  of  gravity  is  an 


TT  ^ 

increase  in  the  distance  of  the  inlet  boundary  with  increasing  F = — . The  positions 

jxU 

of  the  outlet  boundary  are  not  greatly  affected  by  those  of  the  inlet  boundary  except  for 
the  effect  of  gravity  in  the  case  of  F > 0.3.  For  instance,  from  the  known  experimental 

data  used  in  Chapter  3,  f is  approximately  0.02  (when  //o=0.01inch,  and  U=56  ^/sec)- 
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In  practice,  the  effect  of  gravity  has  some  significance  in  the  relationship  between  the 
outlet  boundary  distance,  and  inlet  boundary  distance,  A more  detailed  dynamic 
analysis  would  be  needed  to  yield  the  relationship  between  and  as  affected  by 
boundary  conditions  involving  surface  tension.  However,  as  done  in  the  previous  section 
the  partially  flooded  inlet  boundary  and  the  Swift-Stieber  outlet  boundary  condition  are 
introduced  to  yield  results  for  the  general  features.  The  results  are  obtained  from  the 
Reynolds  equation  with  different  inlet  boundary  distances  and  Swifl-Stieber  outlet 
conditions. 

Equations  (4.29)  and  (4.39)  show  that  the  dimensionless  film  thickness,  0 , only 
depends  on  the  outlet  boundary  distance,  Surface  tension  is  not  accounted  for  in  the 
model  by  Swift-Stieber  (Reynolds)  boundary  condition.  However,  surface  tension  is 
accounted  for  the  model  by  the  separation  boundary  condition  connected  with  Coyne  and 
Elrod  boundary  condition.  Both  cases  in  Figure  (4.3)  uggests  that  the  dimensionless  film 
thickness,  0 , after  rupture  at  the  outlet  zone  should  depend  only  on  the  distance  of  the 
outlet  boundary,  as  shown  in  Figures  (4.4)  and  (4.5).  The  dimensionless  film 

thickness,  0,  as  function  of  the  distance  of  outlet  boundary,  can  be  used  to 
determine  the  distance  of  the  outlet  boundary,  , depending  on  the  Capillary  Number. 

The  film  thickness  after  rupture  show  that  the  case  of  Reynolds  boundary  condition  is 
much  thicker  than  that  of  separation  boundary  condition.  This  relationship  will  provide 
some  help  in  modeling  the  flow  with  high  shear  rate  in  the  lubricating  zone  to  be 
discussed  layer. 

It  is  apparent  that  the  model  based  on  the  lubrication  approximations  is  not 
adequate  possibly  at  small  value  of  due  to  numerical  instability.  The  relationship 

between  ^^and  is  developed  for  explaining  starved  inlet  and  outlet  boundary.  From 
this  phenomenon  and  analytical  solution  the  results  show  the  model  has  limits  on  the  size 
of  the  rotation  system.  This  needs  a reasonable  disc  size  to  give  starved  boundary 
condition  and  to  avoid  numerical  difficulty  in  small  size.  It  seems  likely  that  one  source 
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of  error  arises  from  the  phenomenon  in  the  case  of  having  surface  tension  acting  on  the 

free  surface  where  entry  zone,  and  exit  zone,  4h » ®^ist  around  the  minimum  gap. 

At  higher  speed  and  moderately  loaded  operation  condition  it  is  already 
recognized  that  Swift-Stieber  (Reynolds)  boundary  condition  is  more  suitable  to  use  in 
the  above  case.  Up  to  high  value  of  speed,  the  film  thickness  by  withdrawal  increases  in 
the  free  surface  with  increasing  speed.  Hence,  high  speed  affects  the  surplus  at  inlet.  It 
also  causes  pressure  rise  in  the  fluid  with  a constraint. 

4.7  Summary  and  Conclusion 

A simple  mathematical  formulation  based  on  the  lubrication  model  can  predict  the 
dimensionless  film  thickness  to  be  formed  after  the  outlet  zone  for  both  Swift-Stieber 
boundary  condition  and  separation  boundary  condition.  It  is  a function  of  outlet  zone 
position  and  dimensionless  film  thickness.  The  position  of  the  outlet  boundary  is 
governed  by  the  inlet  boundary  position.  The  position  of  the  inlet  boundary  is  mainly 
determined  by  the  geometrical  parameter  WR,  rotating  speed,  and  Ca. 

The  difficulty  in  the  analysis  lies  in  the  complexity  of  the  flow  in  the 
neighborhood  of  the  regions  of  the  inlet  and  outlet  boundary  zones.  Some  research  points 
out  that  a recirculation  pattern  exists  at  outlet  zone  near  stationary  surface  [16,  17].  It 
would  not  explain  the  conditions  required  for  the  validity  of  a lubrication  model,  in  which 
it  is  assumed  that  the  streamlines  are  nearly  parallel  throughout  the  flow  region. 
Actually,  it  needs  appropriate  boundary  conditions  at  the  separation  point.  The  major 
point  is  that  the  Swift-Stieber  condition  is  not  suitable  to  deal  with  the  existence  of  a 
negative  pressure  before  the  separation  and  recirculating  flow. 
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Figure  4. 1 Cylinder-flat  plate  geometry  and  nomenclature 
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Figure  4.2  Relationship  between  inlet  and  outlet  (^  ) 

position  for  various  values  of  T (gravity  parameter)  by  using 
Swift-Stieber  boundary  condition. 


0=H./H 
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Figure  4.3  Dimensionless  film  thickness  occurred  from  rupture 

region  as  function  of  exit  zone, 

(Comparison  of  Swift-Stieber  B.C.  and  Coyne  & Elrod 
and  Separation  B.C.) 
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Figure  4.4  Relationship  between  inlet  zone  and  outlet  zone  of  the 
lubricant  using  Swift-Stiber  (Reynolds)  boundary  condition 
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Figure  4.5  Relationship  betwen  inlet  zone  and  outlet  zone  of  the 
lubrication  by  using  separation  boundary  condition,  and  Coyne  - 
Elrod  boundary  condition 
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CHAPTER  5 

THERMOHYDODYNAMIC  ANALYSIS  IN  THE  LUBRICATING  REGION 

5.1  Introduction 

When  there  is  oil  present  in  the  vicinity  of  two  surfaces,  a lubricating  film 
manages  to  move  between  the  surfaces  and  completely  separate  them.  This  allows  the 
load  exerted  on  the  surfaces  to  be  supported  or  transmitted  with  small  friction  forces.  The 
lubrication  of  these  surfaces  can  be  achieved  in  different  ways  which  can  be  specially 
shaped  to  ensure  lubrication,  as  in  the  slider  and  cylinder-plane  cases  shown  in  Figure  (1- 
1-b).  This  is  called  hydrodynamic  lubrication  (HDL).  In  this  rolling  or  sliding  element 
systems,  the  combined  effect  of  pressure  and  temperature  on  the  oil  work  together  to 
produce  a film.  This  type  of  lubrication  is  called  thermohydrodynamic  lubrication 
(THD). 

THD  exists  in  successfully  lubricated  non-contacts  between  heavily  loaded  rolling 
and  sliding  machine  components  that  do  not  conform  to  each  other.  This  type  of 
lubrication  is  characterized  by  the  ability  to  generate  pressures  and  temperature  that  are  of 
sufficient  magnitude  to  increase  the  lubricant  viscosity  considerably.  The  combined 
physical  effects  which  constitute  the  complete  THD  problem  are  as  follows: 

1.  The  hydrodynamic  effect  implies  that  the  lubricant  obeys  the  laws  of  viscous 
flow; 

2.  The  thermal  effects  are  significant  when  the  shear  heating  is  large  enough  to 
significantly  alter  the  viscosity  of  lubricant  due  to  the  high  temperature; 
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3.  The  variable  viscosity  effect  is  very  important  especially  in  heavily  loaded  cases 

where  the  pressure  and  temperature  changes  are  substantial. 

The  analysis  of  THD  problem  has  proven  to  be  highly  complex.  One  of  the 
difficulties  arises  from  the  fact  that  the  viscosity  of  the  lubricant  varies  exponentially  with 
temperature.  Another  difficulty  is  setting  its  boundary  conditions  resulting  from 
cavitation  effects.  The  resulting  governing  equations  are  highly  nonlinear  differential 
equations  combining  of  the  hydrodynamic,  and  energy  equations.  The  numerical  solution 
of  these  equations  is  very  difficult  since  its  convergence  is  extremely  sensitive  to  the 
changes  in  the  major  parameters,  especially  the  properties  of  the  lubricant. 

5.2  Literature  Review 

In  the  early  1880s  Beauchamp  Tower  performed  experimental  investigation  to 
study  friction  in  journal  bearings.  Osborne  Reynolds  eventually  introduced  the  concept 
of  the  hydrodynamic  lubrication  theory  and  derived  the  Reynolds  equation  describing  the 
pressure  distribution  in  hydrodynamic  bearings  [22].  The  design  of  hydrodynamic 
bearings  is  generally  based  on  the  assumption  of  constant  viscosity  without  heat  effects. 
Early  works  considering  such  heat  effects  in  fluid  film  system  considered  the  change  of 
the  entire  lubricant  viscosity  with  temperature  rise  in  systems  such  as  bearings.  Cope 
systematically  derived  the  governing  equations  for  viscous  flow  in  bearing  systems  [29]. 
Based  on  the  assumption  that  lubricant  oil  carries  away  the  entire  frictional  heat,  he 
derived  a simplified  energy  equation  in  which  he  neglected  the  variation  of  viscosity 
across  the  film. 

After  Cope’s  work,  which  is  known  as  adiabatic  theory,  some  authors  explained 
the  effect  of  temperature  in  slider  bearings  [26,  27,  30,  32,  48].  They  found  that  the 
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pressure  distribution  based  on  the  adiabatic  theory  is  significantly  lower  than  that 
predicted  by  the  isothermal  theory.  Starting  from  this  theory  Zrenkiewicz  included  the 
temperature  variation  across  the  film  in  the  energy  equation  and  solved  the  governing 
equation  of  a simple  two-dimensional  parallel  surface  bearing  [49].  For  the  variation  of 
viscosity  both  along  and  across  the  fluid  film,  Dowson  first  proposed  the  generalized 
Reynolds  equation  which  is  a function  of  viscosity  based  on  temperature  [50].  A model 
which  predicts  bearing  performance  based  on  realistic  thermal  boundary  conditions  on 
the  stationary  and  moving  surfaces  and  variation  of  the  film  viscosity  with  temperature  is 
finally  introduced  as  the  thermohydrodyanmic  (THD).  With  the  advent  of  high-speed 
digital  computers,  Dowson  and  Hudson  presented  THD  solutions  for  a slider  bearing  in 
which  the  transfer  to  heat  to  the  stationary  solid  was  included  [51-53].  Ezzat  and  Rohde 
also  studied  THD  effects  in  a finite  inclined-plane  slider  bearing  [54-56].  They  solved 
the  Reynolds  equation,  a three-dimensional  energy  equation  for  the  fluid  and  a two- 
dimensional  heat  conduction  equation  for  the  stationary  solid  along  with  a viscosity- 
temperature  relation  for  the  lubricant  simultaneously. 

Using  a similar  numerical  solution  procedure,  Hueber  proposed  a study  of  the 
THD  effects  in  sector  trust  bearings.  His  study  contains  a comparison  of  the  THD  results 
with  isothermal  and  adiabatic  boundary  conditions  [57,  58].  He  concluded  that  bearing 
performance  is  strongly  affected  by  the  inlet  temperature,  speed,  and  film  thickness. 
Most  recently,  Seireg  and  Wang  presented  a complete  THD  study  of  slider  bearings  and 
journal  bearings  in  which  a three-dimensional  variation  of  lubricant  viscosity  and  thermal 
expansion  were  considered  [26,  27].  They  explained  the  occurrence  of  high  temperature 
at  the  central  zone  of  system.  Comparing  their  findings  with  experimental  results,  they 
found  a considerable  difference  leading  to  the  conclusion  that  THD  solution  by  traditional 
shear  zone  may  not  be  completely  satisfactory. 
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5.3  Design  Parameters 

In  the  design  of  rotating  systems  such  as  bearings,  the  responsibility  of  an 
engineer  is  to  provide  customers  with  the  system  geometry,  material  specifications,  and 
method  of  lubricant  supply.  For  a given  design  condition,  a number  of  parameters  are 
engaged  in  determining  the  above  specifications.  To  make  the  design  process  clear,  the 
parameters  can  be  classified  into  some  parts  according  to  their  characteristics,  as  shown 
in  Table  4. 


Table  4:  Design  Parameters  of  Sliding  Systems  (Bearings,  Disc) 


1 

Parameter 

1 Performance 

Friction 

Load 

Lubricant  Flow 

Control 

Maximum  Temperature 
Maximum  Pressure 
Minimum  Film  Thickness 
Speed 

Geometry 

Inclined  Angle  or  eccentricity 
Width  and  Length,  or  Diameter 

Auxiliary 

Lubricant  Supply  Temperature 
Lubricant  Type 

Note:  Newtonian  case 


Variables  which  influence  the  rotating  system  performance  are  as  follows: 

■ Pressure:  Lubricant  pressure  is  developed  due  to  a lubricant  flow  in  the  gap. 
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■ Temperature:  Lubricant  temperature  increases  above  the  supply  temperature  due  to 
viscous  dissipation  of  the  lubricant. 

■ Film  thickness:  The  variation  of  the  film  thickness  determines  the  pressure  and 
temperature  distributions. 

■ Viscosity  of  lubricants:  Viscosity  is  highly  dependent  on  temperature,  and  influences 
the  load  capacity  significantly. 

■ Angle:  slopes  at  which  the  pad  is  declined  with  respect  to  the  runner 

These  variables  are  related  with  one  other,  as  shown  schematically  in  Figure  (5.1). 

5.4  Dependence  of  Viscosity  on  Temperature  and  pressure 

This  section  presents  the  equation  of  state  used  to  relate  temperature  and 
viscosity.  Viscosity  is  the  most  important  property  of  the  oil  or  lubricant  employed  in 
THD,  since  it  is  a determining  factor  in  pressure  and  film  fluid  dynamics  between  the 
rolling  and  sliding  bodies. 

Experiments  have  shown  that  viscosity  is  very  sensitive  to  temperature  and 
pressure  changes[22.  36,  59].  Especially  at  high  sliding  speeds  between  the  surfaces,  both 
pressure  and  temperature  within  the  lubricant  film  rise  significantly  as  the  lubricant 
passes  through  the  lubricating  zone.  Accordingly,  an  approximate  formula  for  the 
dependence  of  viscosity  on  temperature  and  pressure  is  necessary  in  order  to  solve  the 
THD  problem  by  numerical  methods. 

There  are  a number  of  empirical  models  available  for  the  dependence  of  viscosity 
on  temperature.  These  are  all  based  upon  the  observation  that  the  viscosity-temperature 
derivative  depends  primarily  on  the  viscosity.  The  simplest  but  the  least  accurate  model 
is  by  Reynolds  [22]. 
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where 


1 ( 

S = In  — 

Tb~~Ta 


temperature/viscosity  parameter 


B 


To,  Mo  oil  irilot  temperature  and  viscosity,  respectively. 

This  study  employs  this  expression  (5.1)  because  the  effect  of  temperature  is  more 
sensitive  to  viscosity  than  that  of  pressure  in  the  THD  problem  except  for  the 
elastohydrodynamic  lubrication  (EHD)  problem.  It  is  more  convenient  to  use  without 
severe  error  in  this  case.  A more  useful  viscosity  temperature  relationship  was  reported 
in  an  reference  [59] 


(5.2) 


where  the  three  constants  all  have  distinctive  physical  meanings. 

The  viscosity-pressure  derivative  is  also  primarily  dependent  on  viscosity.  Bams 
assumed  a linear  relationship  and  the  following  viscosity  pressure  dependence 


M = Moe‘^  (5.3) 

where  Mo  is  the  viscosity  at  atmospheric  pressure  and  a is  a constant.  This  equation 
becomes  inaccurate  above  0.5  GPa,  thus  making  it  inappropriate  for  this  analysis  [36]. 


5.5  Simplified  Reynolds  Equation 


The  equation  that  governs  the  generation  of  pressure  in  the  lubricant  film  is 
known  as  the  Reynolds  equation.  It  is  the  foundation  for  hydrodynamic  lubrication 
analysis  and  was  first  derived  by  Reynolds  in  1886  [22].  The  Reynolds  equation  can  be 
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derived  by  applying  the  basic  equations  of  motion  and  continuity  to  the  lubricant.  Its 
relevant  form  for  thermal  Newtonian  flow  is  derived  as  follow.  The  generalized 
Reynolds  equation  [50]  which  combines  the  momentum  and  continuity  equations  also 
expresses  the  variation  of  density  and  viscosity  across  and  along  the  lubrication  film  is  as 
follows  : 
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Before  deriving  the  much  simplified  equation,  the  following  simplifying  assumptions  are 


made: 


1 . Body  forces  are  negligible,  due  to  the  negligible  mass  of  the  film 

2.  Lubricant  flow  is  laminar,  since  viscosity  is  high  and  Reynolds  number  is 
small  and  the  lubricant  is  a Newtonian  fluid. 

3.  Velocity  gradients  are  only  significant  across  the  film  thickness,  because  film 
thickness  is  much  smaller  than  the  dimensions  of  the  lubricating  zone. 

4.  Inertia  and  surface  tension  forces  are  negligible  compared  to  viscous  forces, 
since  the  lubricant’s  mass  and  acceleration  are  small. 
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5.  There  is  no  slip  at  the  boundary  surfaces  and  solid  boundaries. 

6.  The  lubrication  of  density  is  considered  to  be  a function  of  temperature  and 
pressure  in  cross  film  direction. 

The  generalized  Reynolds  equation  (5.3.1)  can  be  simplified  as 


where 
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The  appropriate  boundary  conditions  for  Reynolds  equation  Eire 

p{x,0)  = p{x,  L)  = p{p,  y)  = p(b,  y)  = 0 for  slider  bearings 
The  thermal  expansion  term  of  equation  (5.5)  in  the  z direction  by  using  the  equation  of 


state  of  the  fluid  can  be  written  as 
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where  a , and  fi  are  the  thermal  expansion  coefficient  and  bulk  modulus  of  the  lubricant 


respectively  [26]. 


5.6  Energy  Equation 

This  section  presents  the  development  of  a parabolic  energy  equation  for  use  in 
finding  the  temperature  field  in  the  lubricant  film.  The  temperature  distribution  within 
the  lubricant  film  can  be  determined  from  the  solution  of  the  conservation  of  energy 
equation  subjected  to  appropriate  boundary  conditions.  Considering  only  conduction  Eind 
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convection  heat  transfer,  compression  and  shear  work,  and  internal  energy,  the  steady 
flow  energy  equation  for  a eontrol  volume  of  the  lubricant  can  be  written  as 

Q + W,Hear  = " Hi„)  (5.7) 


where  Q is  the  rate  of  heat  flow  into  the  volume,  W shear  is  the  rate  of  shear  work  done  on 
the  volume,  m is  the  mass  rate  of  flow  through  the  volume,  and  H is  the  enthalpy  of  the 
fluid. 


Before  preceding  further,  the  following  simplifying  assumptions  are  made; 

1.  Velocity  gradients  are  significant  only  across  the  film  thickness  (as  in  the 
development  of  the  generalized  Reynolds  equation) 

2.  The  density  of  lubricant,  p , within  the  THD  film  is  constant. 

3.  The  film  is  mueh  shorter  in  the  z-direction  compared  to  x and  y-directions  so 
that  the  gradient  of  the  veloeity  component  in  the  z direetion  is  ignored. 
However,  temperature  gradients  are  not  ignored  in  those  directions. 

4.  The  thermal  eonductivity  of  the  lubricant,  k , within  the  THD  film  is  constant. 

5 . Conductions  in  x,  y and  z-directions  are  not  neglected  eompared  to  eonduction 
to  boundary  surfaces. 

Using  the  above  assumptions,  the  energy  equation  for  incompressible,  thin  film  used  in 
this  study  is 
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The  form  of  the  energy  equation  (5.8)  represents  heat  generated  both  due  to  shear  and 
compression  when  a fluid  elements  is  subjected  to  a veloeity  difference  across  its  top  and 
bottom  surfaces.  This  heat  is  dissipated  into  the  solid  surfaces  whieh  the  film  contacts  by 
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conduction,  and  convection.  The  energy  equation  is  subject  to  the  initial  condition  as 
follows; 


T[x,y,z;0)  =Tt„ 


(5.9) 


and  the  appropriate  boundary  conditions  for  the  energy  equation  are 
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= 0 : adiabatic 


= 0 : THD 


(5.10) 


5.7  Dimensionless  Analysis  for  the  Governing  Equations 


There  are  many  advantages  to  non-dimensionalizing  the  governing  equations. 
One  is  that  the  results  have  generality  and  the  problem  of  using  different  system  of  xmits 
is  removed.  Another  big  advantage  is  that  the  number  of  external  variables  is  reduced, 
from  those  having  dimensions  to  their  various  combinations  in  dimensionless  group. 

The  following  dimensionless  variables  are  used  to  nondimensionalize  the 
governing  equations  [26] : 
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Using  the  above  variables  the  Reynolds  equation  can  be  written  in  dimensionless  form  as 
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The  dimensionless  boundary  conditions  for  Reynolds  equation  are 


p{x\0)  = p 
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(5.13) 


The  dimensionless  form  of  the  energy  equation  is  [26] 
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where  [54] 
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The  dimensionless  initial  conditions  for  the  energy  equation  is 


T'[x\y  ,z\o)  = 1 


(5.15) 


and  the  appropriate  non-dimensional  THD  boundary  conditions  for  the  energy  equation 
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The  fluid  film  velocity  profiles  can  be  obtained  from  momentum  equations  and  expressed 


as  a function  of  pressure,  and  viscosity  as  follows: 
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Boundary  conditions  for  Equation  (5.17)  are 


u = U,v  = 0 eX  z = 0 
u = 0,v  = 0 at  z = h 


The  velocity  profile  can  be  expressed  as: 
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Using  the  transformation  variables  (5.1 1),  the  dimensionless  velocity  profiles  are  given 
by 
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The  cross  film  velocity  w may  be  obtained  from  the  continuity  equation  and  integration 
with  w = 0 at  z = 0 and  h as  follows 
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5.8  Performance  Parameter  of  the  Lubricating  Systems 


The  following  parameters  play  a major  role  on  the  hydrodynamic  performance. 
The  load  carrying  capacity  for  the  slider  moving  system  is; 


LB 
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(5.23) 


The  friction  force  can  be  computed  from, 
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The  inlet  flow,  outlet  flow,  and  side  leakage  are. 
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There  is  no  side  leakage  flow  in  the  cavitation  region  x > Xcmuauon  since  the  axial  pressure 

gradient,  — , is  zero.  Due  to  cavitation,  the  film  is  incomplete  and  the  system  is  not 
dz 


covered  with  lubricant  over  the  length. 
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Figue  5.1  Relationship  between  variables  in  the  lubrication  problem  of  the  sliding  motion 


CHAPTER  6 

NUMERICAL  TECHNIQUE 


6.1  Introduction 

Finite  element  methods  (FEM)  have  been  widely  used  in  analysis  of  several 
research  fields  because  irregular  geometry  shape  can  be  easily  accommodated  at  the 
element  level  without  any  difficulty  for  the  setup  of  programming.  Also,  the  size  of  the 
element  can  easily  be  varied  over  the  region  allowing  the  use  of  small  elements  whem 
necessary.  This  method  has  been  gradually  applied  to  deal  with  the  problem  of  film 
thickness  variation  in  lubrication  analysis  such  as  elastohydrodynamic  (EHL)  problem. 

The  FDM  where  the  physical  domain  of  a lubricant  film  can  effectively  be 
transformed  to  a rectangular  computational  domain  is  selected  due  to  its  suitability  to 
THD  problems  without  severe  film  thickness  variation  [26,  27,  60-62].  The  equation  of 
motion  for  a homogeneous  fluid  in  the  absence  of  a finite  rate  chemical  reaction  or  mass 
diffusion  is  based  on  three  physical  conservation  laws.  Those  are  conservation  of  mass 
(continuity),  conservation  of  linear  momentum  (Newton’s  second  law),  and  conservation 
of  energy  (first  law  of  thermodynamics).  First,  the  Reynolds  or  energy  equation  is 
decoupled  fi’om  the  system  of  equations  composed  of  the  continuity  and  momentum 
equations.  Therefore,  the  velocity  and  pressure  fields  are  computed  initially,  and 
subsequently  the  energy  equation  may  be  solved  for  the  temperature  distribution. 
Second,  due  to  the  reduction  of  unknowns  within  the  domain,  the  computer  storage  and 
memory  requirements  are  reduced.  The  numerical  solution  of  the  THD  problem  requires 


76 


77 


the  simultaneous  solution  of  the  Reynolds,  film  shape  (especially  important  in  the  EHL 
problem),  load,  energy,  eind  viscosity  equations.  The  Reynolds  and  energy  equations  are 
formed  to  deal  with  constituent  elements  using  the  central  finite  difference  method. 
These  approximated  equations  are  then  solved  by  numerical  methods.  For  this  class  of 
problems  where  the  governing  equations  can  be  solved  sequentially,  it  is  advantageous  to 
write  a finite  difference  formulation  suited  for  a staggered  grid. 

6.2  Finite  Difference  Form  of  the  Reynolds  Equation 

For  numerical  computation,  a finite  difference  form  of  the  Reynolds  equation  is 
introduced.  The  dimensionless  Reynolds  equation  is 
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Using  the  following  expressions  (6.2),  and  the  central  differences  approximations  and 
substituting  the  Equation  (6. 1 ) 
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Hence,  equation  (6.1)  can  now  be  written  in  finite  difference  form  with  the  above 
substitutions  as: 
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In  equation  (6.3),  O, , O2 » O3 , and  04  refer  to  the  coefficients  applied  to  the  nodes 
- around  the  central  node,  O5 , and  Oe  is  the  source  term. 

The  following  two-step  method  is  used.  Equation  (6.4)  is  solved  iteratively  using 
the  point  successive-over-relaxation  (SOR)  scheme  with  optimal  relaxation  factor,  co . 
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where  p is  the  intermediate  calculated  pressure  value  and  k is  the  index  for  iteration.  In 
order  to  stop  the  output  loop  iteration,  the  following  convergence  criterion  is  used  [26] : 
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where  ERR  is  varied  between  0.02  and  0.03,  depending  on  the  required  accuracy  of  the 
solution. 


6.3  Finite  difference  form  of  the  Energy  Equation 


The  numerical  scheme,  called  alternating  direction  implicit,  is  employed  for 
solving  this  parabolic  energy  equation.  The  ADI  belongs  to  a class  of  methods  known  as 
“Approximate  Factorization.”  In  these  methods,  the  original  multidimensional  FDEs  are 
replaced  by  a series  of  finite  difference  equations  with  can  be  represented  by  tridiagonal 
formulations.  The  advantage  of  the  ADI  scheme  lies  in  its  simplicity  and  efficiency.  The 
tridiagonal  matrix  inversion  is  very  efficient  to  carry  out  on  any  computer.  The  implicit 
scheme  computation  removes  the  stringent  time  step  limit.  The  dimensionless  energy 
equation  (5.18)  can  be  expressed  as  [26] 
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By  center  difference  scheme  and  the  Taylor  series,  equation  (6.6)  becomes  as  follows: 
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Assuming  / = = /" , ADI  scheme  is  as  follows: 
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Applying  the  proper  difference  formula  and  making  element  equation  (6.1 1),  the  resulting 
equations  (6.12),  (6.13)  can  be  solved  with  Thomas  algorithm. 
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In  terms  of  the  format  to  take  advantage  of  the  tridiagonal  form  of  the  coefficient  matrix, 
Equations  (6.12)  and  (6.13)  can  be  rewritten  as 


bT  ]\i-]  +dT  N +aT  i^+ 1 = c 


(6.14) 


where  a,  b,  and  c are  coefficients. 
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>t>  * 

For  boundary  conditions,  Tg  is  known  at  one  boundary  and  r at  the  other.  In 


the  range  of  1 to  N-1,  all  known  T s are  collected  in  to  the  c term.  The  system  can  be 


cast  into  the  following  form  although  the  first  and  last  equations  in  the  array  may  result 
from  the  auxiliary  relationship. 
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Referring  to  the  tridiagonal  matrix  of  coefficients  above,  the  system  is  put  into  an  upper 
triangular  form  by  computing 
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then  computing  the  unknowns  from  back  substitution  and  then 
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6.4  Grid  System 

Figures  (6.1)  and  (6.2)  illustrate  the  finite  moving  system  configuration  and 
geometry.  Diagrams  of  journal  and  slider  bearings  are  sketched  with  nomenclature.  A 
journal  bearing  consists  in  its  simplest  form  of  two  rigid  cylinders.  The  outer  cylinder  is 
usually  held  stationary,  while  the  inner  cylinder  rotates  at  an  angular  velocity.  A slider 
bearing  consists  of  two  inclined  plane  surfaces.  Film  thickness  in  the  lubricating  region 

can  be  established  on  the  angular  coordinate,  0,  for  a journal  bearing  and  on  the  x 
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coordinate  for  a slider  bearing.  By  the  dimensionless  transformation,  a mesh  of 
rectangular  cells  with  constant  sizes  for  the  full  shear  zone  model  and  variable  sizes  for 
the  partial  shear  zone  model  is  established  for  the  development  of  the  finite  difference 
approximation.  The  lubricating  region  consists  of  the  inlet,  and  outlet  as  shown  in  Figure 
(6.1)  for  a journal  bearing  and  Figure  (6.2)  for  a slider  bearing.  This  is  shown  in  Figure 

(6.3).  in  which  5x  , 3y  , and  dz  express  the  element  size  for  x , y , and  z 
respectively.  The  number  of  grids  could  be  fixed  for  different  operating  conditions.  The 
dimensions  across  the  film  thickness  and  film  flow  length  are  divided  into  60  grids.  Film 
flow  width  is  divided  into  10  grids.  The  structure  of  the  finite  difference  mesh  is  shown 
on  Figure  (6.4).  The  geometrical  configuration  of  generalized  lubricating  zone  illustrated 
in  Figure  (6.4)  is  expressed  as  two  structures  with  thin  shear  zone  and  full  shear  zone. 
The  computational  domain  of  finite  moving  system  is  plotted  in  two  ways,  one  having  a 
full  shear  zone  for  the  lubricating  region,  the  other  having  three  different  shear  zone  for 
stationary,  and  moving  regions  and  a central  thin  shear  zone  to  be  discussed  in  Chapter  7. 

A variable  mesh  nodal  structure  like  the  one  shown  in  Figure  (6.3-b)  is  used  to 
provide  close  spacing  in  and  around  the  thin  shear  zone  to  be  discussed  in  next  chapter. 
In  this  figure,  the  inlet  is  to  the  left  and  the  outlet  is  to  the  right.  The  variable  mesh  helps 
to  minimize  the  errors  that  can  occur  because  of  the  large  gradients  in  the  thin  shear  zone 
with  high  shear  rate. 


6.5  Consistency  and  Stability 

As  seen  in  Equation  (6.11),  the  numerical  scheme  for  each  of  the  two  fractional 
steps  from  the  energy  equation  is  relatively  simple.  However,  it  is  necessary  to  show  that 
the  combination  of  these  two  steps  is  consistent  with  the  original  governing  equation. 
This  is  very  hard  to  prove  theoretically  because  of  the  existence  of  the  three  dimensions 
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and  two  additional  heating  source  terms.  This  is  a weak-point  in  applying  the  method  of 
fractional  steps  to  the  problems  of  fluid  dynamics  in  several  variables.  However,  this 
weak-point  does  not  appear  to  affect  the  accuracy  of  the  numerical  model.  For  numerical 

accuracy,  it  is  desirable  to  keep  the  time  step  sufficiently  small,  (typically  At  =2.5x10'^ 

sec)  such  that  Courant  number  ( |v|  = |c  At/ Ax|  < 1 ),  where  c is  coefficient)  does  not 

become  much  lager  than  1 . 

6.6  Computational  Procedure 

The  computational  procedures  are  summarized  as  follows.  The  flow  chart  for  the 
computer  program  is  shown  in  Figures  (6.5).  At  the  start  of  the  numerical  simulation  the 
pressure  and  temperature  fields  are  set  to  the  inlet  lubricant  pressure  and  temperature. 
The  velocity  distribution  of  the  film  is  calculated  at  each  time  step  with  the  continuity  and 
Newtonian  shear  velocity  gradient  relationship  as  shown  in  Figure  (6.5).  However,  when 
the  shear  zone  ratio,  h^/h,  is  evaluated  as  a function  of  operating  system  characteristic 
number  to  be  introduced  in  the  next  chapter,  the  program  requires  proper  initial  velocity 
values  for  each  zone  of  this  model  as  shown  in  Figure  (6.3).  Also,  it  structures  all  shear 
action  takes  place  in  the  shear  zone  which  is  at  the  center  of  the  film. 

(1)  Supply  all  the  physical  input  data.  The  data  are  the  rotating  system  geometry,  oil 

« 

t 

type,  reference  oil  temperature,  operating  speed,  and  properties  of  the  lubricants. 

Values  of  the  input  data  allow  the  determination  of  nondimensional  parameters. 

(2)  Initialize  pressure,  velocity,  and  temperature. 

(3)  Generate  uniform  grids  of  x,  y,  z coordinates. 

(4)  Compute  coefficients  for  the  simplified  generalized  Reynolds  equation  and  solve 

the  simplified  generalized  Reynolds  equation  to  obtain  a new  pressure  field  for  the 


film. 
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(5)  Check  the  pressure  profile.  If  the  convergence  criterion  is  not  satisfied,  then  the 
current  pressure  profile  is  adjusted  by  iterative  computation  until  it  is  satisfied. 

(6)  Calculate  the  velocity  profile  in  the  full  shear  zone. 

(7)  Update  viscosity  and  determine  temperature  profile  using  energy  equation  with  its 
boundary  conditions. 

(8)  Repeat  loop  (4)  to  (8)  until  the  pressure  and  temperature  distributions  reach  steady 
state. 

(9)  Compute  the  load  from  the  pressure  distribution,  friction  force,  flow  rate,  and  heat 
convection.  Output  the  pressure,  the  temperature  distribution,  and  the  bearing 
performance  data. 


6.7  Numerical  Calculation  and  Comparison 

An  extensive  amount  of  the  thermohydrodynamic  results  are  presented  for  slider 
bearing  systems.  Slider  bearing  systems  in  typical  hydrodynamic  cases  are  chosen  for 
analyzing  thermohydrodynamic  phenomenon.  These  solutions  are  obtained  using  the 
finite  difference  scheme  described  in  the  previous  chapter.  A computer  code  is  written  in 
Fortran  language  to  obtain  numerical  solutions  to  the  THD  problem.  This  code  is 
originally  developed  by  Wang  to  determine  the  shear  zone  ratio  which  are  compatible 
with  the  experimental  data  performed  for  journal  bearings  [62].  For  the  results  presented 
in  the  present  study,  his  algorithm  is  mainly  applied  to  give  a more  realistic  solution  to 
the  problem. 

The  program  is  capable  of  generating  results  for  any  bearing  geometry,  and 
lubricant  properties.  The  important  design  parameter  and  test  fluid  for  the  test  slider 
bearings  are  presented  in  Table  5 for  which  other  numerical  results  are  available.  The 
solutions  are  validated  by  comparing  them  to  published  experimental  data  and  numerical 
results.  This  chapter  then  proceeds  to  develop  a generalized  process  for  THD  analysis. 
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First,  the  validity  of  the  energy  equation  with  compression  heating  source  is  examined  by 
comparing  the  steady  state  distribution  in  temperature  and  pressure  to  those  published  in 
the  literature.  Based  on  the  validated  algorithm,  a new  method  is  developed  to  use  a 
rheological  concept  model  to  the  adiabatic  shear  band  and  the  glass  transition  behavior 
for  comparison  with  recent  experimental  results. 


Table  5:  Design  Parameters  for  Test  Slider  Bearings  [26,  27,  62] 


Case 

# 

r Length, 
B (inch) 

Width, 
L (inch) 

^min 

(inch) 

(inch) 

Velocity 

(ft/sec) 

SAE 
(oil  #) 

Temp 

(”F) 

1 * 

3 

3 

0.002 

0.004 

100 

30 

100 

2 * 

3 

3 

0.002 

0.004 

50 

30 

100 

3 * 

3 

3 

0.002 

0.004 

200 

30 

100 

4 * 

3 

3 

0.002 

0.004 

100 

50 

100 

< * 

3 

3 

0.002 

0.004 

100 

30 

150 

6 * 

3 

3 

0.002 

0.004 

100 

10 

100 

Note:  *:  Full  shear  zone  model  (Newtonian  case) 


Table  6:  Lubricant  viscosity-temperature  data,  ')  [26, 62] 


SAE  oil  # 

5 

10 

30 

50 

ju , lb- sec/ ft^ 

3.24  X 10-" 

6.6384  X lO"* 

0.00221 

0.00504 

S 

-0.0142 

-0.0172 

-0.022 

-0.025 
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Table  7:  Common  Input  Parameters  on  Bearing  [26,  62] 


1 k = 0.075 Btu  / {hr  - ft 

-°F)  : conductivity 

p = 54.6lb„,l  ft^ 

: density 

Cp  = 0.Amullb„,-F 

: specific  heat 

a = 0.00036"  r‘  : 

thermal  expansion  coefficient 

yff  = 31.68x10^6// ft^ 

: bulk  modulus 

6.8  Model  with  Full  Film  Shear  Zone  in  the  Lubrication  Zone 

Developments  in  lubrication  technology  have  been  widely  dedicated  to  industrial 
areas.  The  journal  bearing  is  still  the  most  encountered  in  engineering  applications. 
However,  for  simple  checks  at  different  speeds  and  wide  range  of  lubricants,  the  slider 
bearing  with  different  geometries  is  considered. 

The  survey  of  such  geometries  is  given  in  the  Table  5.  One  of  the  main 
consideration  is  that  viscous  dissipation  in  the  lubricating  film,  separating  two  surfaces  in 
relative  motion  can  results  in  high  rates  of  heat  generation.  The  purpose  of  this 
comparison  is  to  confirm  the  numerical  methods  in  the  Newtonian  case  by  comparing 
them  to  the  numerical  results  of  Ezzat  and  Rohde  [54].  The  system  used  for  analysis  is 
the  three  dimensional  slider  bearing  systems  shown  in  Figure  (6.2).  The  length  to  width 
ratio  is  1 .0.  An  isothermal  and  adiabatic  analysis  are  considered  for  the  stationary  solid 
surface.  Transition  motion  of  the  slider  bearing  is  calculated  and  is  compared  with  Ez2at 
and  Rohde’s  results  [55].  The  pressure  and  temperature  traces  of  the  slider  bearing  center 
are  plotted  in  Figure  (6.6)  throughout  Figure  (6.9).  Based  on  the  steady  state  solution 
after  500  time  steps  (0.01  sec),  the  pressure-time  and  temperature-time  traces  for  a 
constant  speed,  100  ft/sec,  SAE  30  oil,  and  100  °F  initial  temperature  are  shown  in  Figure 
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(6.6)  as  typical  time  history  results.  The  fluid  film  pressure  is  shown  to  increase 
proportionally  during  the  first  few  milliseconds  of  operation,  almost  reaching  its  steady- 
state  value.  This  is  the  transient  time  to  reach  constant  running  speed  from  motionless 
situation.  Next,  the  transient  response  of  the  fluid  film  under  various  operating 
conditions  is  investigated.  Results  are  presented  graphically  in  Figure  (6.7)  throughout 
Figure  (6.9).  Another  pressure-time  and  temperature-time  plots  for  a constant  speed  (100 
ft/sec)  is  shown  in  Figure  (6.7)  for  different  lubricant  oils.  Both  pressure  and  temperature 
values  increase  with  time  and  are  relatively  high  for  higher  number  lubricant  oil  since  the 
initial  value  of  viscosity  associated  with  viscous  heat  generation  is  higher  under  the 
conditions.  The  time  scale  is  of  the  order  of  milliseconds.  The  effect  indicates  the  higher 
the  running  speed,  the  shorter  the  time  to  reach  steady  state. 

Figures  (6.8)  and  (6.9)  represent  the  results  after  applying  isothermal  condition  to 
the  stationary  pad.  The  effects  of  the  initial  temperature  is  presented  in  Figure  (6.8).  The 
effects  of  running  speed  on  the  transient  performance  is  also  shown  in  Figure  (6.9).  As 
shown  in  Figure  (6.11),  the  higher  maximum  temperature  corresponding  to  the  adiabatic 
boundary  condition  are  responsible  foe  the  decreased  maximum  pressure  compared  to  the 
isothermal  boundary  condition.  Pressure  values  in  the  hydrodynamic  lubrication  regime 
with  the  very  thin  high  shear  effects  can  be  obtained  by  keeping  the  generated  heat  in  the 
central  thin  lubricating  zone.  This  will  be  discussed  in  the  next  section.  The  rheological 
results  with  central  thin  adiabatic  heat  band  and  high  shear  rate  are  also  discussed.  Both 
steady-state  solutions  of  maximum  pressure  and  maximum  temperature  in  the  transient 
performance  are  in  close  agreement  with  the  numerical  solution  by  Ezzat  and  Rohde  [54, 
55].  The  main  reason  there  for  the  trivial  difference  which  exists  between  the  two  results 
is  due  to  the  three-dimensional  analysis  used  and  to  the  addition  of  the  compression 
heating  source  in  the  energy  equation  in  the  present  study. 
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The  pressure  distribution  in  the  fluid  film  is  probably  the  most  important  factor  in 
the  analysis  of  the  performance  of  lubricating  system.  The  basic  characteristics  of  all  the 
pressure  profiles  are  a thermohydrodynamic  pressure  built  up  in  the  lubricating  zone  and 
a pressure  drop  to  the  atmospheric  pressure  in  the  inlet  and  outlet  zones.  Figure  (6.10) 
shows  a typical  pressure  distribution  in  the  lubricating  zone  of  60  x 10  x 60  mesh  grid 
with  side  leakage.  The  pressure  profiles  in  the  lubricating  zone  for  various  operating 
conditions  (running  speed)  are  shown  in  Figure  (6.11).  The  pressure  distribution  is 
strongly  influenced  by  the  running  speed  and  the  thermal  boundary  condition  in  the  fluid 
film.  The  location  of  maximum  pressure  is  slightly  changed  when  using  different  thermal 
boundary  conditions. 

Figures  (6.12-b),  (6.13),  and  (6.14)  are  the  corresponding  isotherm  contours  in  the 
central  plane  (defined  as  the  y=L/2  plane).  The  effect  of  the  thermal  boundary  condition 
at  the  stationary  pad  is  shown  in  Figure  (6.12).  These  results  show  very  similar  trends  as 
the  results  of  Ezzat  and  Rohde.  In  the  case  of  the  adiabatic  condition  no  heat  is 
transmitted  to  the  stationary  pad  and  a relatively  high  temperature  boundary  layer  is 
formed  as  shown  in  Figure  (6.12-a).  Using  both  surface  temperatures  of  the  bearing  as 
the  initial  temperature  has  the  effect  of  reducing  the  temperature  in  the  fluid  film  as 
shown  in  Figure  (6.12-b).  Most  of  heat  dissipated  is  conducted  away  to  the  stationary 
pad  causing  the  maximum  temperature.  Maximum  temperature  is  only  160°F  in  Figure 
(6.12-b).  Otherwise,  Figure  (6.12-a)  shows  the  temperature  contours  when  the  interface 
between  stationary  pad  and  adjacent  lubricant  is  assumed  to  be  adiabatic.  From  the 
contours,  it  can  be  seen  that  there  is  small  temperature  variations  along  the  film 
lubrication  except  around  inlet  zone.  This  boundary  condition  is  unrealistic.  But,  it  does 
serve  to  emphasize  the  influence  of  the  thermal  boundary  condition.  It  also  indicates  the 
probability  of  the  existence  of  such  high  temperature  adiabatic  shear  boundary  layer  in 
the  central  thin  lubricated  region.  Figures  (6.13)  and  (6.14)  show  the  effect  of  running 


90 


speed  and  lubricant  oil  type  on  the  fluid  film  thermal  transients.  In  all  cases,  the 
stationary  pad  is  assumed  to  be  isothermal. 

Figure  (6.15)  shows  the  velocity  profiles  in  a full  shear  zone  in  the  direction  of 
motion.  In  the  velocity  stream  the  effect  of  the  thermal  boundary  and  running  speed  is 
demonstrated  by  the  different  shear  rates  at  the  region  where  the  maximum  pressure  is 
acting.  This  steady-state  THD  performance  for  inclined  slider  bearing  under  the 
assumption  of  a full  shear  zone  in  fluid  film  is  performed  without  the  effect  of  heat 
conduction  to  the  stationary  pad  in  fluid  film. 

6.9  Summary  and  Conclusion 

This  chapter  presents  a detailed  numerical  method  that  is  employed  for  the 
prediction  of  THD  performance  of  fluid  film  bearings.  A two-dimensional  Reynolds 
equation  for  the  pressure  distribution,  three-dimensional  full  energy  equation  for  the 
temperature  distribution,  and  velocity  equation  are  transformed.  The  governing  equation 
which  combines  the  three  main  equations  is  solved  in  two  ways  for  full  shear  and  central 
shear  respectively.  The  regular  grid  system  is  used  for  the  full  shear  zone  model.  Some 
of  the  results  using  this  approach  are  used  to  validate  the  developed  3 dimensional 
analysis  with  the  assumption  of  lubrication  with  full  shear  zone  by  comparing  the 
obtained  results  with  published  numerical  results. 
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Figure  6. 1 Diagram  sketch  of  journal  bearing 
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Figure  6.2  Diagram  sketch  of  a slider  bearing 
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Figure  6.3  Nodal  structure  of  computational  domain  of 

lubricating  zone 
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Figure  6.4  Finite  difference  mesh 
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Figure  6.5  Program  flow  chart  for  THD  analysis  using  the  Newtonian  behavior  and 

full  shear  zone  in  lubricating  region 
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Figure  6.6.  Slider  bearing  pressure  and  temperature 
history  for  full  shear  zone  (SAE  30  oil,  Tj^=100  °F, 

U=100  ft/sec,  B=L=3.0  inch,  h|^jj^=0.002  in,  hj^g^=0.004  in) 
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Figixre  6.7  THD  pressure  (a),  temperature  (b)  - time  history 
Tjj^  =100  °F,  U=100  fit/sec,  slider  bearing,  B=L=3.0  in, 
hjjjjjj=0.002  in,  h^g^=0.004  in 
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Figure  6.8  THD  temperature  -Time  history 
SAE  30  oil,  U=100  ft/sec 

Slider  bearing  (B=L=3.0  in,  hj^jjj=0.002  in,  hj^^j^=0.004  in) 


Maximum  Temperature, 
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Figure  6.9  THD  temperature-Time  history 
SAE  30  oil,  Ti„=100  °F 

[Ezzat  and  Rohde]  - steady  state  THD  temperature  solution 
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Figure  6.10.  Pressure  distribution 

U=100  ft/sec,  Tj„=100  ®F,  SAE  30  oil,  B=L=3.0  in, 

hmin=0.002  in,  h^g^=0.004  in 
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Figure  6.1 1 Pressure  distribution  along  the  middle  plane  of 
the  slider  bearing  in  the  direction  of  motion 

(SAE  30  oil,  Tjjj=100  °F,  B=L=3  in,  h„,i„=0.002  in,  h„ax=0-004  in) 


axis  TT  - axis 


Figure  6.12  Mid-plane  dimensionless  film  temperature  distribution 
U=100  ft/sec,  Tjj^^lOO^F, 

B=L=3.0  in,  SAE  30  oil,  hjjjjjj=0.002  in,  h^a^=0.004  in) 

(a)  adiabatic  boundaries  for  upper  stationary  pad 

(b)  boundaries  with  isothermal  for  upper  pad 


axis  TT  - axis 
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Figure  6.13.  Dimensionless  film  temperature  distribution  of 

the  bearing  in  the  direction  motion 

(B=L=3.0  in  slider  bearing,  hj^jj^=0.002  in,  hj^^g^=0.004  in, 

Tin=100  ®F,  U=100  ft/sec,  ) 

(a) SAE  10  oil 

(b)  SAE  50  oil 


axis  77  - axis 
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Figure  6.14  Dimensionless  film  temperature  distribution  of  the 
slider  bearing  in  the  direction  of  motion 

B-L=3.0  in,  h^j„=0.002  in,  hj^g^=0.004  in,  SAE  30  oil,  Tj„=100  °F, 
fixed  temperature  boundary  condition  for  upper  surface  pad 


(a)  U=50  ft/sec 

(b)  U=200  fl/sec 


axis  TT  - axis 
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Figure  6.15  Mid-plane  dimensionless  film  velocity  distribution 
U=100  ft/sec,  Tjj^=100  ^F,  B=L=3.0  in  slider  bearing,  SAE  30 

(a)  adiabatic  boundaries 

(b)  boundaries  with  isothermal  Tjj^  for  stationary  pad 


CHAPTER  7 

THE  RHEOLOGICAL  MODEL 


7.1  Introduction 


The  shear  characterization  of  the  thin  high-pressure  lubricant  film  has  been  the 
subject  of  investigation  in  recent  years  [26,  27,  63,  64].  The  shear  rheological  behavior 
of  the  liquid  under  severe  operating  conditions  displays  several  unique  and  interesting 
features  that  shear  theory  predicts  and  the  liquid  may  undergo  a glass  transition.  Starting 
with  the  work  of  Milne  many  works  have  been  performed  to  evaluate  the  effect  that  non- 
Newtonian  lubricants  may  have  on  the  performance  of  bearings  [65-71].  Some  works 
listed  on  Table  8 show  the  diversity  of  lubricant  models  employed  and  the  different 
effects  produced  by  the  various  analytical  approaches  [72].  Such  analytical  approaches 
and  experimental  data  can  reduce  the  uncertainties  prevailing  in  this  area.  Some  studies 
in  EHD  condition  have  shown  that  the  lubricant  behavior  caimot  be  described  by  a 
classical  Newtonian  model  [64,  66,  67,  69-71,  73].  Also,  some  experiments  support  the 
unique  rheological  behavior  which  is  observed  in  high  pressure  EHD  lubrication 
conditions  [64,  74].  Some  models  for  describing  the  rheological  behavior  and  thermal 
effects  in  lubricated  EHD  contacts  have  been  proposed  to  fit  the  experimental  finding 
[64,  70,  71].  Also,  shear  bands  have  been  also  observed  in  a high  temperature  flow 
visualization  tests  [64,  74]. 

In  general,  the  development  of  rheological  models  considering  the  isothermal 
shear  deformation  in  a homogeneous  liquid  is  not  likely  to  result  in  significant  advances 
in  the  field.  In  lubricating  systems  with  EHD  contacts,  relatively  simple  methods  have 
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been  proposed  based  on  experimental  results  [73].  Wang  and  Seireg  have  proposed  a 
method  to  simulate  the  shear  behavior  of  the  film  [27].  They  simulated  the  rheological 
effects  by  classifying  three  different  shear  zones  in  Newtonian  fluid  films  [26].  Their 
results  showed  the  significant  thermal  effects  in  the  central  thin  shear  zone  and  indirectly 
verified  the  existence  of  the  thin  shear  zone  through  correlations  with  experimental 
results.  However,  their  model  cannot  describe  all  experiments  and  control  some  effects 
of  the  rheological  trend  because  of  using  uniform  shear  zone  ratio  along  the  lubrication 
flow. 

In  this  chapter  the  special  model  that  is  used  in  the  thermohydrodynamic  analysis 
of  lubrication  problem  is  introduced.  It  incorporates  a modified  velocity  equation 
characterized  by  the  inclusion  of  the  non-Newtonian  rheological  model.  It  will  also  be 
based  on  the  Bingham  plastic  model  which  is  one  of  the  rheological  models  describing 
the  shear  rheological  response  of  the  lubricant  in  the  lubricating  region.  Theoretical 
studies  of  lubrication  with  Bingham  model  have  explained  the  existence  of  liquid  fluid  (a 
possible  rigid  part)  to  be  attached  to  one  or  the  other  surface  by  authors  [75,  76].  A 
model  of  the  behavior  of  Bingham  fluids  is  presented  by  Tichy  to  exhibit  a yield  stress 
[77,  78].  He  used  a continuous  mathematical  function  as  the  constitutive  law  and  a 
generalized  Reynolds  equation  to  solve  general  lubrication  problems  and  also  indicated 
the  influence  of  the  stress  parameter.  This  study  will  deal  with  viscous  phenomena  such 
as  the  case  of  a stationary  layer  and  a uniform  moving  layer  by  introducing  the  concept  of 
the  yield  stress.  Furthermore,  the  goals  of  this  study  are  (i)  to  present  a continuous 
formulation  which  is  applicable  through  the  flow  field  to  characterize  the  behavior  (ii)  to 
calculate  the  behavior  in  such  lubrication  conditions,  and  (iii)  to  compare  these  results 
with  the  predictions  of  the  empirical  theory  (Square  root  relationship  between  pressure 
and  velocity  in  a thin  film)  where  applicable. 
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Table  8:  Some  Works  on  Non-Newtonian  Lubricants  [72] 


Year 

Author 

Lubricant  Model 

Bearing  to  Application 

1957 

Milne 

Viscoelastic 

Slider 

1967 

Hsu 

Pseudoplastic 

Journal 

1973 

Tipei  & Rohde 

Directional  Viscosity 

Slider 

1978 

Hamoy 

Viscoelastic 

Journal 

1983 

Bourgin  & Gray 

Polymer  Thickened 

Journal 

1984 

Buckholz 

Rheological 

Journal 

1991* 

Tichy 

Bingham  Plastic 

Journal 

1992* 

Tichy 

Bingham  Plastic 

Slider 

Note:  * refer  references  [77,  78]. 


7.2  Non-Newtonian  Behavior 


A generalized  form  of  Newton’s  law  can  be  expressed  as: 

T = riy  (7.1) 

where  the  viscosity  t]  is  any  function  of  the  scalar  invariants  of  the  shear  rate  f . If  the 
non-Newtonian  viscosity,  rj , is  dependent  on  the  components  of  the  tensor  f then  it  also 
depends  on  combinations  of  the  components  of  the  scalar  tensor. 

The  three  independent  combinations  are 


7i  = Sr//,  h = TTfyfjn  and  h = yf  jkf u 

i i j i j k 


(7.2) 


The  first  invariant  of  the  strain  rate  vanishes  for  an  incompressible  fluid, 
/,  = 2(v-r)=0.  For  shearing  flows,  the  third  invariant  of  strain  rate  becomes  zero. 
Hence,  the  viscosity,  rj,  will  he  considered  to  depend  only  on  I2  [79-81].  For  the 
lubrication  geometry  illustrated  in  Figure  (2.1),  w « v and  djdy  « djdz , the  shear  rate 
magnitude,  7 , can  be  expressed  as  follows: 
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(7.3) 
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(7.4) 


. du  dv 

where  v^,  = — and  = — . 

dz  ^ dz 

The  components  of  the  shear  stress  along  the  direction  of  film  flow  and  film  width  flow 
are; 
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(7.5) 


Form  Equation  (7.5)  the  non-Newtonian  viscosity  for  Bingham  plastic  model  can  be 
defined  as; 


Y 

Therefore,  the  stress  components  can  be  expressed  as  follows; 


(7.6) 


.\du 


Txz  = rjif)-r 
dz 


dz 


(7.7) 


du  dv 

In  most  lubrication  problems,  — » — . Hence  Equation  (7.7)  becomes  as  follows 

dz  dz 


Txz  = M^  + u and  Xyz  = viy)^^^ 
dz  dz 


(7.8) 


Figure  (7.1)  shows  some  examples  of  the  shear  rate  and  yield  stress.  In  the  first 
group  the  fluids  are  without  yield  stress,  ro  > curve  (a)  is  for  a Newtonian  fluid  for  which 
the  viscosity,  rj,  is  independent  of  the  shear  rate,  f • Curve  (b)  describing  the  shear 
thinning  fluid  shows  increasing  shear  stress  less  than  in  proportion  to  the  shear  rate,  f • 
Curve  (c)  describing  the  shear  thickening  fluid  shows  increasing  shear  stress  in 


proportion  to  the  shear  rate,  y • la  the  second  group  of  fluids  with  yield  stress,  xo  > the 
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fluid  flows  when  the  shear  stress  exceeds  a yield  stress.  Otherwise,  the  fluid  behaves  as 
an  elastic  solid  when  the  shear  stress  does  not  exceed  the  yield  stress.  Curve  (d) 
describes  Bingham  plastic  model.  Curves  (e)  and  (f)  are  for  fluids  behaving  non-linearly 
with  a yield  stress.  The  Bingham  plastic  model  of  non-Newtonian  fluid  behavior  has 
been  used  to  describe  the  behavior  of  lubricants  such  as  greases,  the  flow  of  melts  and 
slurries.  Theoretical  studies  of  lubrication  with  the  Bingham  model  will  be  based  on  the 
concept  of  rigid  layers  adhering  to  the  surfaces  at  all  locations  along  the  fluid  film. 

7.3  Rheological  Model  Description 

Some  lubricants  exhibit  a critical  shear  behavior  which  can  be  modeled  as  a 
Bingham  plastic  under  severe  condition.  The  ideal  Bingham  plastic  model  is 
characterized  by  two  parameters,  a yield  stress  and  a viscosity  dependent  on  stress  [78]. 
In  this  model,  the  fluid  behaves  as  a rigid  solid  when  the  magnitude  of  the  stress  is  less 
than  the  yield  stress.  Rigid  layers  can  occur  in  the  fluid  film  as  three  different  types, 
stationary,  uniform  moving,  and  shear  according  the  pressure  gradient  in  the  sliding 
direction. 

In  order  to  understand  the  non-Newtonian  effect  of  lubricants  by  applying  the 
concept  of  Bingham  plastic  model  into  the  THD  problem  with  constant  film  thickness  it 
is  necessary  to  develop  an  appropriate  Reynolds  type  equation  and  modified  elliptical 
type  velocity  equation  which  can  handle  the  formation  of  layers  in  the  fluid  film.  When 
stationary  and  uniform  moving  layers  exist  in  the  fluid  film.  Equation  (2. 1 9)  representing 
the  velocity  is  no  longer  valid.  Therefore,  a generalized  Reynolds  equation  and  an 
appropriate  rheological  model  must  be  considered.  Among  the  governing  equations 
general  descriptions  for  the  velocity  equation  with  rheological  results  are  introduced.  It  is 
necessary  to  explain  the  reasons  for  employing  the  best  elliptical  shaped  velocity  profile 


110 


of  rheological  behavior  in  the  lubrication  flow  between  two  surfaces.  Studies  on  non- 
Newtonian  fluids  have  mainly  depended  on  the  models  that  is  employed  to  decide  on  a 
proper  choice  the  constitutive  equation.  This  means  totally  different  type  governing 
equations  for  the  problem  under  these  conditions,  specially  velocity  profile.  As  the 
differential  and  rate  type  have  none  of  these  conditions  for  the  non-Newtonian 
phenomenon  of  lubrication  problem,  none  of  these  models  can  adequately  describe  the 
experimental  results  that  have  been  established. 

The  starting  point  in  the  flow  characteristics  of  the  non-Newtonian  fluid  under 
consideration  from  the  flow  of  a Newtonian  fluid  is  the  existence  of  narrow  layer  around 
central  lubrication  region  where  the  velocity  gradients  are  exceedingly  high.  The 
existence  of  this  thin  layer  is  known  to  take  place  in  certain  non-Newtonian  fluids  [64]. 
This  phenomenon  can  not  be  explained  in  terms  of  the  classical  method  of  the  differential 
type  from  the  Newtonian  case.  To  avoid  this  difficulty  the  present  model  characterizes 
two  different  and  antisymmetric  profile  shapes  which  velocity  gradients  are  sharply 
changed  in  the  central  zone  and  the  high  pressure  region  to  describe  the  rheological 
properties  of  non-Newtonian  fluids.  Conclusively,  the  main  assumptions  on  which  the 
velocity  equation  is  based  to  explain  the  rheological  properties  at  the  high  shear  rates  and 
shear  stress  at  the  interface  between  the  two  layers  and  the  shear  zone  can  be  listed  as 
follows: 

1 . Shear  stress  is  proportional  to  the  velocity  gradient  across  the  film. 

2.  The  fluid  can  exhibit  non-Newtonian  behavior  under  certain  extreme 
lubrication  condition  (thin  films,  high  pressure  etc). 

3.  The  fluid  will  not  slip  at  the  interfaces  between  the  two  layers  and  the  shear 
zone. 

4.  The  pressure  across  the  film  is  assumed  to  be  constant. 

5.  The  film  thickness  is  very  small  and  thus  the  flow  is  assumed  to  be  laminar. 
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6.  At  the  center  location  across  fluid  film  the  velocity  gradient  is  assumed  to 
change. 

7.  ta  and  Tc  in  Figure  (7.2)  have  low  shear  stress  values  along  the  fluid  film  in 
the  full  shear  zone.  And,  td  and  te  have  the  minimum  shear  stress  values 
for  the  yield  stress  along  the  film  flow. 

8.  Two  oils  of  the  same  viscosity  at  a given  test  temperature  and  pressure  are 
considered  for  the  upper  region  and  lower  region  in  the  shear  zone. 

9.  It  is  assumed  that  velocity  in  the  central  region  is  a half  of  the  uniform  moving 
velocity. 

10.  The  lubricant  is  assumed  to  be  a Newtonian  fluid. 

1 1 . For  the  lower  region  in  the  shear  zone,  the  Newtoni2m  velocity  equation  of  the 
differential  type  is  employed  with  the  proper  boundary  condition  including 
assumption  (9).  While,  for  the  upper  region  in  shear  zone,  velocity  profile  has 
the  anti-symmetric  shape  of  the  lower  region. 


When  the  layers  form  within  the  lubricant  film  zone  they  are  assumed  to  move  at 
velocities,  Uuppet  = 0.  Ucemmi  = and  ui„wer  = U along  the  direction  of  motion.  The 

velocity  profiles  and  their  corresponding  shear  stress  in  the  xy  plane  can  be  illustrated  as 
shown  in  Figure  (7.2).  Point  A and  C indicate  a shear  stress  ta  and  tc<to-  The  point 
on  the  centerline  indicates  a shear  stress  which  is  larger  than  a yield  stress.  Referring  to 
Figures  (7.2),  and  (7.3)  the  boundary  conditions  for  the  velocities  can  be  written  as: 
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. = U V „ = 0,M  ,=0,v  =0 

z=0  z=0  ’ z=h  ’ z=0 
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= U , u =0,  u 
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' = /jc= 


^=h, 


= 0 


(7.9) 


In  the  lower  region  of  the  central  shear  zone,  z < ho  the  velocity  distributions  and 
their  gradients  are; 
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(7.10) 

(7.11) 


du 

dz 


(7.12) 
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dz  rj  dy\_  Fl 


(7.13) 


he  ^ he  2^  he^  he  ^ 

where  Fo  = I ~dz , F\  = \ —dz , Fl  = j—dz , and  F\  = j—dz  . 

h,n  HlV  0 7 0 7 

The  cross-film  velocity  w can  be  obtained  by  using  the  continuity  Equation  (2.14)  and 

two  differential  forms  with  respect  to  x,  and  y,  respectively,  for  the  lower  shear  region. 
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(7.14) 


In  the  upper  region  of  central  shear  zone,  hc^z<h„,  the  velocity  distributions  and  their 
gradients  are; 


1 ,,  It/  , \ U A ^ dp 
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(7.16) 

(7.17) 

(7.18) 
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The  velocity  distributions  and  their  gradients  in  the  rigid  layer  of  0<z<fiL  are, 


u = U , V = 


du 

dz 


and  in  the  rigid  layer  region  of  are, 


(7.20) 


o dudv 
u = 0,  v = — = — = 0, 

dz  dz 


(7.21) 


A Reynolds  equation  for  this  type  of  fluid  having  layers  (stationary  layer  and 


uniform  moving  layer)  can  be  derived  by  using  Equation  (7.10)  throughout  Equation 
(7.21)  into  the  continuity  equation  and  integrating  the  resulting  equation  across  the  fluid 
film.  The  result  is  given  in  the  similar  form  of  Equation  (5.5)  obtained  in  Chapter  5.  A 


generalized  Reynolds  equation  that  employs  this  model  is  given  below. 
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where 
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A similar  set  of  equations  as  those  of  the  previous  Newtonian  case  are  developed 
for  the  pressure  behavior  of  the  lubricant  under  severe  operating  condition  with  a 
Newtonian  motor  oil.  In  the  case  of  no  solid  layer  attached  to  upper  stationary  surface 


Fa 


encountered  in  the  Bingham  model,  the  parameters  F*  and  — f disappear. 


Fo 


To  account  for  the  property  variation  with  temperature,  the  energy  equation  must 
be  considered.  A form  of  energy  equation  for  the  non-Newtonian  fluid  can  be  written  as 
follows. 


Pcp-—  = '^  •{kVT)-¥f^v  (7.23) 

With  proper  formulation  of  the  boundary  conditions,  a three-dimensional  energy  equation 
taking  into  account  the  temperature  variation  in  the  sliding  direction,  side  direction,  and 
across  the  lubricant  film  yields  remarkably  accurate  results  for  THD  problems. 

In  the  adiabatic  band  having  very  high  temperature  due  to  high  shear  rate  at  the 

central  zone,  0 < < //^ , the  simplified  energy  equation  and  boundary  conditions 

associated  with  the  energy  equation  are; 
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(7.24) 


The  boundary  conditions  associated  with  the  energy  equation  for  the  shear  zone  are 


expected  as; 


T>,^,y,,z)=To 
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(7.25) 


In  the  uniform  velocity  moving  layer,  Q<z\‘^h\,  the  energy  equation  becomes; 


^z\ 


+pr, 

\u^] 

dx\_ 

(7.26) 


and  the  boundary  conditions  for  the  above  region  are, 

ri(0,>^i,zi)=  ri(xi,y|  ,0)=7’o 


(7.27) 


In  the  stationary  layer  zone,  Q<  zi^  hi,  the  energy  equation  and  associated 


boundary  conditions  are. 
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7.4  Non-dimensionalization  of  the  Governing  Equations  for  Rheological  Behavior 


For  the  model,  the  governing  equations  and  boundary  conditions  are  non- 
dimensionalized  using  the  following  parameters  in  addition  to  the  parameters  used  in  the 
Newtonian  case  as  follows; 
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(7.30) 


The  non-dimensional  viscosity  for  the  present  model  as  in  the  Bingham  plastic  model  can 
be  written  as; 


(7.31) 
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The  dimensionless  form  of  Equation  (7.7)  which  is  valid  at  the  central  shear  zone, 

(O  < z,  < hX  is; 
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Otherwise,  at  the  other  zones. 
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(7.34) 


In  the  lower  region  0<Zi  ^^/^l  of  the  flow  in  the  shear  zone,  the  non- 


dimensional  velocity  distribution  and  gradients  are; 
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the  computational  domain  of  the  shear  zone. 
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In  the  upper  region 


z.  ^ h 


of  the  flow  in  the  shear  zone,  the  non-dimensional 


velocity  distribution  and  gradients  are 
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The  non-dimensional  velocity  distributions  and  their  gradients  in  the  rigid  layer  region  of 
0<zi<hi  are, 


_ , _ 5m  dv 
M =1,  V = — = — = 0, 

dz\  dz\ 


(7.45) 


and  in  the  rigid  layer  region  of  0 < ^2  ^ /j  2 are, 


_ ^ _ 5m  5v 
M = 0,  V = — = — = 0, 

dzi  5z2 


(7.46) 


A non-dimensional  Reynolds  equation  for  this  type  of  layered  fluid  film  (stationary  layer, 
uniform  moving  layer,  and  shear  zone)  can  be  similarly  applied  to  Equation  (6.1)  of 


Chapter  6 for  the  two  dimensional  case  in  the  x-y  plane. 
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In  the  adiabatic  band  having  very  high  temperature  due  to  the  high  shear  rate  at 
the  central  zone,  0<Zs-hs^  the  simplified  non-dimensional  energy  equation  is; 
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The  non-dimensional  boundary  conditions  associated  with  the  energy  equation  for  the 
shear  zone  can  be  considered  as; 


(7.49) 


In  the  uniform  moving  velocity  layer,  0 < ^ /j  i , the  energy  equation  becomes; 
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and  the  boundary  conditions  for  the  above  region  are, 

fi(0.5^pzi)  = f,(jci,j,,0)  = l 


(7.50) 
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In  the  stationary  layer  zone,  0 < Z2^  hi’  the  non-dimensional  energy  equation 


and  associated  boundary  conditions  are. 
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and  the  boundary  conditions  for  the  stationary  region  are. 
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7.5  Numerical  Scheme 

The  dimensionless  equations  are  discretized  and  solved  by  the  finite  difference 

method.  A forward  formula  is  used  for  — while  central  difference  approximation  is 

dt 

used  for  all  other  derivatives.  The  Reynolds  equation  is  an  elliptic  partial  differential 
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equation  and  is  solved  using  the  Successive  Over  Relaxation  method  and  the  energy 
equation  is  solved  by  a fully  implicit  scheme.  The  difference  form  of  the  Reynolds 
equation  for  the  present  model  is  already  introduced  in  Chapter  6.  The  Successive  Over 
Relaxation  method  of  two  steps  can  be  used  in  the  same  manner  because  an  improved 
estimate  of  the  pressure  values  is  computed  by  applying  the  formula  of  Equation  (6.14)  at 
every  grid  point  without  considering  the  effect  of  the  non  constant  grid  size  in  the  across 
film  direction. 

In  the  shear  zone,  0<zs^hs>  the  non-dimensional  energy  equation  can  be 
expressed  as  the  following  element  equation; 
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where 

A,=^,  — + A2  = ^2  — + <^2^^and  = ^ — + ,f,  — 

ox  dx  oy  dy  oz  dz 

In  the  uniform  moving  zone,  0 < z,  < ^, , the  ADI  scheme  for  the  non-dimensional  energy 
equation  is  as  follows; 
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In  the  uniform  stationary  zone,  0 < Z2  ^ 2 ’ the  ADI  scheme  for  the  non-dimensional 
energy  equation  is  as  follows; 
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7.6  The  Computational  Procedure 

The  shear  rheological  behavior  of  the  lubricant  in  THD  lubrication  has  been  the 
subject  of  considerable  research.  The  lubricant  under  highly  transient  conditions  is 
subjected  to  relatively  severe  operating  conditions.  High  pressure  and  high  temperature 
at  the  central  zone  can  induce  a complex  rheological  response  of  lubricant.  In  the  liquid 
region  the  shear  rheoregical  behavior  would  be  expected  to  be  a classical  viscous 
behavior  with  possible  rheological  behavior  at  high  rates  of  change  of  stress.  All  the 
dimensionless  governing  equations  and  boundary  conditions  are  discretized  using  the 
finite  difference  method.  To  predict  the  non-Newtonian  behavior  of  the  present  model 
based  on  the  Bingham  plastic  model,  the  most  important  step  in  addition  to  the  previous 
Newtonian  case  is  to  determine  the  boundaries  of  the  uniform  moving  layer,  stationary 
layer,  and  shear  layer  in  the  lubricating  zone.  Figure  (7.4)  shows  the  program  flow  chart 
for  THD  calculation  of  the  non-Newtonian  fluid  films. 

An  arbitrary  velocity  profile  where  the  velocity  gradient  is  changed  at  the  central 
line  of  the  film  thickness  is  assumed  in  the  entire  shear  zone  as  in  the  Newtonian  case  in 
Chapter  6 for  the  initial  condition.  The  corresponding  pressure  and  temperature  profiles 

are  also  initiated  properly.  In  this  analysis,  no  rigid  layer  formation,  hi  = 0 in  the 

uniform  moving  velocity  zone  and  = ^ in  the  stationary  velocity  zone,  is  assumed  for 
the  first  computational  cycle.  The  corresponding  magnitude  of  the  shear  rate  is  computed 
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from  the  initial  velocity,  pressure,  and  temperature  profiles.  Next,  the  Reynolds  equation 
is  solved  using  the  SOR  method.  Having  computed  the  pressure  distribution,  the  velocity 
is  computed  and  then  the  energy  equation  for  the  fluid  is  solved.  At  every  cycle,  the 
pressure,  temperature,  velocity,  and  shear  zone  ratio  are  updated  until  the  steady  state  is 

reached.  Here,  the  boundaries  of  h\  and  hi  in  the  film  thickness  are  determined  from 
computed  shear  stresses  and  the  resultant  shear  stress  at  every  cycle.  These  points  along 
the  lubrication  film  are  determined  to  define  the  boundaries  of  the  regions  where  the 
resultant  shear  stress  is  less  than  the  yield  stress.  Calculated  zones  are  updated  by 
Newzone=01dzone-01dzone*l/(NX*5)*2  regions.  For  a typical  lubricating  system 
simulation,  the  number  of  mesh  points  in  the  axial  direction  is  set  to  Ay  = 1 0 , A;,  = 60 

mesh  points  along  the  lubrication  film,  and  Az  = 60  ( A,  = 20 , A = 20 , and  A = 20 ) 

*“  Zs  Z\  Z2  ' 

mesh  points  across  the  film  thickness. 

The  following  steps  are  the  procedure  to  determine  the  existence  of  zones  having 
uniform  moving  and  stationary  layers  in  the  lubricating  zone. 

1 . Assume  that  there  is  full  shear  formation  in  the  lubricating  zone  for  the  first 
computational  cycle. 

2.  Compute  the  resultant  shear  stress  with  initial  velocity  distribution. 

3.  Solve  the  Reynolds  equation  in  the  full  shear  zone  assumed  in  step  1 . 

4.  Calculate  the  velocity  distribution  of  the  non-Newtonian  behavior  by  using  the 
model  with  the  gradient  changing  at  the  central  zone  across  the  film  thickness. 

5.  Compute  the  resultant  shear  stress  and  compare  it  with  the  yield  stress  at 
every  nodal  points  ( Az  = 60 ) across  lubrication  film  thickness  and  centrally 
axial  points.  Determine  the  region  where  the  resultant  shear  stress  is  less  than 
the  yield  stress.  Once,  shear  zone  ratio  is  computed,  the  values  are  updated  to 
proceed  next  step. 
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6.  Having  determined  the  new  shear  zone  ratios,  step  2 ~ step  5 is  repeated  until 
the  computational  cycle  reaches  the  final  cycle  for  the  steady  state. 

7.  Print  out  all  performance  characteristics  of  the  film. 

7.7  Empirical  Evaluation  of  the  Thin  High  Shear  Zone  and  Experimental  Data 

The  basic  concept  of  shear  zone  model  is  that  shearing  is  assumed  to  occur  only 
in  a small  zone  located  centrally  in  the  fluid  film  instead  of  making  the  assumption  of 
shear  in  the  entire  film  [26, 27].  Due  to  the  high  rate  of  heat  generation  in  the  shear  zone, 
the  film  is  thought  to  expand  in  the  across  film  direction  contributing  to  the  pressure  rise 
in  the  film.  The  key  to  this  concept  lies  in  the  time  scale  of  the  fluid  expansion  relative  to 
the  time  scale  of  the  fluid  flow  through  the  bearing  space.  The  basic  idea  is  that  the  fluid 
expansion  rate  is  much  faster  than  the  flow  rate  which  results  in  pressure  development 
[31]. 

A constant  shear  zone  thickness  ratio  was  assumed  and  the  model  was  solved  for 
a fixed  bearing  geometry  and  running  conditions.  This  was  repeated  until  the  peak 
pressure  determined  by  the  model  matched  the  experimental  peak  pressure  for  the  UW-1 
bearing  for  the  same  operating  conditions.  Upon  checking  the  pressure  distributions 
predicted  by  the  model,  both  axially  and  circumferentially,  it  was  found  that  the 
predictions  were  nearly  exact.  This  appeared  to  confirm  the  original  hypothesis  of  a 
shear  zone.  This  process  was  repeated  for  several  experimental  conditions  from  which  an 
empirical  relationship  predicting  the  ratio  of  the  shear  zone  height  the  film  height  was 
determined.  The  ratio,  hs/h  was  found  to  be  a function  of  the  product  Sommerfeld 

number  S and  the  Peclet  number,  Pe,  — = f{S-  p^)  [27]. 
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where  Sommerfeld  and  Peclet  numbers  are  defined  as  follows. 
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where  p is  the  average  pressure  at  the  film  and  N is  the  speed  in  the  journal  bearing. 

The  behavior  of  an  isoviscous  fluid  film  wedge  can  be  characterized  based  on 
Newtonian  fluid  dynamics  by  Sommerfeld  number.  The  transient  behavior  of  heat 
transfer  in  the  fluid  wedge  is  characterized  by  the  Peclet  number.  To  develop  the 
empirical  shear  zone  relationship  the  UW-1  bearing  was  selected  due  to  the  extensive  test 
data  available  for  various  conditions  [62].  The  test  conditions  for  the  present  model  was 
selected  from  experiments  which  cover  speeds  from  100  to  2400  rpm  for  journal  bearings 
and  from  10  to  40  ft/sec  for  slider  bearings.  The  lubrication  oil  ranged  from  SAE  10  to 
50,  and  the  eccentricity  ratio  for  journal  bearing  ranged  from  0.6  to  0.9  with  inlet 
temperatures  from  100  to  160  °F.  Using  this  empirical  relationship  for  h/h  the  model 
was  checked  by  available  experimental  data.  In  all  cases  the  model  predictions  matched 


well  with  the  experimental  data.  The  results  from  the  model  predicted  the  maximum 
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pressures  and  the  pressure  distribution  in  all  cases.  To  sum  up  empirical  evaluation  of  the 
shear  zone,  the  test  conditions  used  are  given  in  Tables  9 and  10. 


Table  9:  Design  Parameters  for  Investigated  Journal  Bearings  (UW-1)  [26, 27,  62] 


Case 

# 

Diameter, 

D(in) 

Width, 

L(in) 

C,  (in) 

8 

Speed 

(rpm) 

SAE 
(oil  #) 

Temp  1 
(°F) 

A 

3.154 

2.0 

0.0134 

0.9 

1050 

10 

128 

B 

3.154 

2.0 

0.0134 

0.9 

2100 

10 

128 

c 

3.154 

2.0 

0.0134 

0.87 

1200 

50 

160 

D 

3.154 

2.0 

0.0134 

0.87 

2400 

50 

160 

E 

3.154 

2.0 

0.0134 

0.6 

1000 

30 

107 

F 

3.154 

2.0 

0.0134 

0.6 

2000 

30 

107 

Table  10:  Design  Parameters  for  Investigated  Slider  Bearings  (UW-3)  [26,  27,  62] 


Case 

Length, 

B(in) 

Width, 

L(in) 

iimin 

(in) 

hmax 

(in) 

Velocity 

(fVsec) 

SAE 
(oil  #) 

Temp  1 
(°F) 

G 

3 

3 

0.004 

0.008 

40 

10 

100 

H 

3 

3 

0.004 

0.008 

20 

10 

100 

I 

3 

3 

0.004 

0.008 

10 

30 

100 

J 

3 

3 

0.004 

0.004 

10 

30 

100 

K 

3 

3 

0.006 

0.006 

10 

30 

100 

L 

3 

3 

0.008 

0.008 

10 

30 

100 
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7.8  A Parametric  Investigation  of  the  Effect  of  the  Rheological  Properties 
7.8.1  Journal  Bearing 

The  results  of  all  cases  of  journal  bearing  (UW-1)  listed  in  Table  9 are  shown  in 
Figure  (7.5)  throughout  Figure  (7.18).  Figure  (7.5)  and  (7.6)  show  the  change  of 
maximum  pressure  with  respect  to  the  rotating  speed  of  the  UW-1  journal  bearing.  For  a 
bearing  speed  at  a specific  operation  condition,  the  shear  yield  stress  parameter  used  was 
obtained  from  the  comparison  with  experimental  data  at  the  high  speed  (N>2000  rpm). 
The  results  of  these  THD  analysis  with  the  thermal  effects  of  the  thermal  expansion  in  the 
across  film  direction  show  the  same  trend  with  the  observed  empirical  relationship 
between  the  operating  speed  and  film  pressure  of  the  bearing.  The  pressure  in  a fixed 
geometry  fluid  film  is  found  to  be  proportional  to  the  square  root  of  speed  [30,  32].  It  is 
interesting  to  note  that  the  journal  bearing  pressure  and  speed  characteristics  of  two 
different  oil  types  show  the  effect  of  lubricant  properties  such  as  viscosity,  inlet 
temperature  and  the  lubricant  film  thickness  as  defined  by  the  eccentricity  ratio  on 
journal  bearing  performance  by  comparing  Figures  (7.5)  and  (7.6). 

Figures  (7.7)  and  (7.8)  show  the  pressure-time  history.  The  results  were  used  to 
iteratively  evaluate  the  yield  shear  stress  parameter  in  the  steady  state  solution  which  best 
fits  the  experimental  pressure.  The  computer  code  is  designed  to  evaluate  the  shear  zone 
across  the  film  flow  as  determined  by  the  shear  yield  stress  parameter  after  a transient 
time  for  reaching  the  operating  bearing  speed.  From  such  transient  time,  the  results  show 
that  the  maximum  pressures  in  the  lubricating  region  are  greatly  influenced  by  the  shear 
yield  stress  parameters.  This  shear  yield  stress  parameter  is  also  an  important  factor  in 
determining  the  shear  zone  ratio  which  leads  to  significant  increase  in  temperature  in  the 
central  narrow  shear  zone.  In  Figure  (7.7),  when  the  bearing  speed  is  2400  rpm,  the 
transient  time  to  reach  this  speed  is  0.025  sec  and  the  steady  state  pressure  is  reached 
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approximately  after  0.5  sec.  The  time  step  used  in  this  numerical  simulation  is  0.001  sec 
for  computational  stability.  The  maximum  pressure  at  the  steady  state  was  used  to 
determine  the  value  of  shear  yield  stress  parameter  which  corresponds  to  given 
experimental  data.  Also,  for  some  given  operating  conditions,  several  shear  yield  stress 
values  were  attempted  in  same  manner.  The  best  shear  yield  stress  parameter  was 
applied  to  find  the  shear  zone  ratio,  hs/h,  by  checking  the  yield  region  of  the  resultant 
shear  stress  along  the  film  flow.  The  varying  shear  zone  ratios  for  the  investigated  case  B 
are  plotted  in  Figure  (7.9).  Figures  (7.10),  (7.1 1)  and  (7.12)  show  the  trend  of  the  shear 
zone  ratio  along  the  film  flow  for  case  B,  D,  and  F of  the  UW-1  journal  bearing.  From 
the  computed  results  in  each  case,  the  dimensionless  shear  zone  along  the  film  flow  are 
found  to  be  in  the  same  order  of  magnitude  as  the  empirical  boundary  shear  layer 
obtained  by  Equation  (7.58)  based  on  the  isoviscous  theory.  A value  of  the  shear  stress 

parameter,  fo=2.8,  is  needed  for  reaching  yield  in  the  solid  boundaries,  A eind  B in 

Figure  (7.2).  This  result  with  fo<2.8  is  also  the  same  as  that  obtained  by  the  THD 
analysis  with  full  shear  zone  in  this  case.  The  calculated  shear  zone  ratio  mainly  depends 
on  the  resultant  shear  stress  which  is  affected  by  the  velocity  gradient  and  viscosity  by 

temperature.  The  value  of  f o = to  hoj Mo  ^ is  mainly  dependent  on  the  value  ro  and  U. 

Also  because  the  objective  of  this  problem  is  to  determine  the  thin  shear  zone  in  the 
central  region  it  can  be  assumed  that  the  other  regions  have  the  characteristics  of  solids, 

were  zo  = P'lA.  As  shown  in  Figure  (7.9)  increasing  the  value  of  shear  yield  stress 
parameter  tends  to  cause  the  formation  of  a wide  region  in  which  the  velocity  is  nearly 
constant,  which  would  be  non-shearing  fluid  in  the  narrow  lubrication  zone.  Under  anti- 
symmetric velocity  behavior  in  the  shear  zone,  the  points  along  the  film  flow  direction 
can  be  systematically  obtained  by  checking  where  the  resultant  shear  stress  at  the  solid 
boundary  is  less  than  the  shear  yield  stress  parameter. 
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The  velocity  profile  of  lubricant  film  at  various  locations  along  film  flow 
direction  are  shown  in  Figure  (7.13)  for  different  operating  conditions  with  SAE  50,  Tin 
=160  °F,  and  eccentricity  0.9.  In  the  beginning  zone  of  lubrication  region  the  velocity 
profile  is  almost  linear.  As  the  gap  narrows  in  the  converging  region  the  lubricant  gets 
squeezed  in  the  thin  shear  zone  and  is  forced  to  change  its  velocity  so  that  the  velocity 
profiles  become  concave  in  each  half  region  of  shear  zone.  The  results  came  from  the 
considered  velocity  equation  which  is  assumed  to  cause  large  shear  stress  at  the  central 
region  in  the  shear  zone  due  to  the  high  shear  rate  generated  by  the  velocity  gradient. 

In  Figures  (7.14)  and  (7.15),  the  pressure  distribution  in  the  fluid  film  show  some 
of  the  most  important  bearing  performance  characteristics  such  as  the  load  carrying 
capacity,  and  the  frictional  force.  The  pressure  distribution  is  strongly  influenced  by 
thermal  effects  in  the  fluid  film.  The  thermal  effects  also  depend  on  the  central  thin  shear 
zone  which  leads  to  significant  temperature  rise.  Figure  (7.14)  shows  that  the  thin 
adiabatic  shear  zone  resulting  from  the  large  shear  yield  stress  parameter  significantly 
reduces  the  pressure  in  the  lubricating  film.  Figure  (7.15)  also  shows  that  the  pressure 
distribution  along  the  center  line  of  journal  bearing  is  significantly  influenced  by  bearing 
speed.  As  shown  in  Figures  (7.14)  and  (7.15),  the  speed  and  the  shear  yield  stress 
parameter  are  the  main  factors  in  determining  the  pressure  distribution. 

The  temperature  distribution  is  displayed  in  the  central  sheeir  zone  and  at 
stationary  and  moving  layers  under  the  yield  zone  with  the  assumed  rheological  velocity 
field.  In  the  thin  shear  zone,  the  zone  is  considered  as  an  adiabatic  shear  band  causing 
high  temperature  in  the  lubricating  region  Both  boundaries  of  the  shear  zone  are 
adiabatic.  For  the  stationary  layer,  the  interface  with  the  stationary  solid  is  assumed  to  be 
insulated.  For  the  moving  layer,  the  temperature  at  the  interface  with  the  moving  solid  is 
assumed  to  have  the  same  value  as  the  lubricant  inlet  temperature  and  the  solid-like  fluid 
is  also  assumed  to  be  insulated.  Figures  (7.16),  (7.17)  and  (7.18)  show  temperature 
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distribution  in  the  fluid  film  at  the  mid-plane  of  the  journal  bearing.  The  shear  zones 
shown  in  the  above  Figures  are  magnified  to  have  the  same  dimension  as  the  other  two 
layers  for  detailed  temperature  profile  in  illustration  of  the  actual  thin  thickness.  Figure 
(7.16)  shows  the  results  at  the  conditions  of  SAE  10  oil  and  high  eccentricity  ratio,  s=0.9. 
It  should  be  note  that  most  of  the  heat  generation  occurs  in  the  central  shear  zone  with 
high  shear  rate.  Around  20“  and  180°,  temperature  profile  show  relatively  high 
temperature  compared  to  any  other  circumferential  location  in  the  journal  bearing.  It  is 
considered  to  be  due  to  the  low  values  of  the  shear  zone  ratio.  Especially,  the  smallest 
value  of  shear  zone  ratio  around  1 80°  in  the  circumferential  location  shows  the  maximum 
temperature  as  can  be  seen  in  Figure  (7.12).  In  the  fluid  film,  of  the  bearing  systems,  the 
moving  layer  has  low  temperature  but  high  gradient  across  film  flow.  There  is  some  heat 
convection  to  the  moving  layer  zone  because  lubricant  in  the  moving  layer  is  quickly 
replaced  by  the  fresh  lubricant  attached  to  the  moving  surface. 

7.8.2  Slider  Bearings  (Tapered  and  Petrallel  Bearings! 

The  results  of  the  G,  H,  I cases  for  the  tapered  slider  bearing  (UW-3)  listed  in 
Table  10  are  shown  in  Figure  (7.19)  throughout  Figure  (7.33).  Also,  the  results  of  the  J, 
K,  L cases  for  the  parallel  slider  are  plotted  in  Figures  (7.34)  to  (7.39).  Figure  (7.19) 
shows  that  the  slider  bearing  pressure  versus  speed  characteristics  also  has  same  trend  as 
the  pressure  in  a fixed  geometry  fluid  film  which  is  proportional  to  the  square  root  of 
speed.  Figures  (7.20),  (7.21)  and  (7.34)  show  the  film  pressure  and  temperature  history 
of  slider  bearings.  From  the  results,  it  can  be  seen  that  the  speed  of  bearing  is  the  main 
key  to  influence  the  pressure  and  temperature  generated  in  the  lubricating  region  of  the 
slider  bearing  of  a fixed  geometry.  As  shown  in  Figure  (7.34),  the  film  thickness 
becomes  also  another  important  factor  in  a parallel  slider  bearing.  The  average  film 
thickness  of  the  tapered  slider  bearing  (hmin=0.004  inch,  hmax=0.008  inch)  is  set  to  be 
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equal  to  the  film  thickness  of  the  parallel  surface  slider  bearing  which  is  h=0.006  inch. 
Both  sliders  have  a square  stationary  surface  3.0  inch  in  length  and  running  speed  at  10 
ft/sec.  The  lubricant  is  SAE  30  and  the  inlet  temperature  is  100  °F.  As  shown  in  Figure 
(7.36),  the  results  of  the  pressure  distribution  along  the  film  flow  direction  shows 
significant  difference  between  cases  I and  K.  It  is  considered  as  the  effect  of  thermal 
edge  in  the  converging  area  which  influences  the  build  up  of  much  higher  pressure.  The 
effect  of  the  thermal  expansion  is  very  insignificant  in  building  pressure  in  the  lubricating 
region  after  reaching  the  operating  speed  of  bearing  as  shown  in  the  pressure  history  of 
Figure  (7.21).  The  yield  parameters  which  determine  shear  zone  ratio  along  film  flow 
direction  has  considerable  influence  on  the  pressure  distribution,  as  shown  in  Figures 
(7.22),  (7.23),  and  (7.35).  The  effect  of  various  yield  parameters  in  the  cases  of  G,  H, 
and  K is  illustrated  by  the  trend  of  shear  zone  ratio  in  Figures  (7.24),  (7.25),  (7.27), 
(7.29),  and  (7.37).  Figures  (7.25)  and  (7.27)  show  the  trend  of  yield  parameters  and  shear 
zones  in  the  process  of  seeking  the  best  agreement  in  the  maximum  pressure  with 
experimental  data  and  other  THD  analysis.  Figures  (7.26),  (7.28)  and  (7.30)  show  the 
thickness  of  the  dimensionless  shear  zone  existing  at  central  region  across  film  flow 
against  the  dimensionless  film  thickness.  The  shear  zone  ratio  for  various  positions  along 
the  film  flow  are  plotted  along  with  the  empirical  shear  zone  ratio  by  Equation  (7.58).  At 
the  value  shear  stress  parameter  determined  for  the  best  correlation  with  the  experimental 
data  the  computed  shear  zone  ratios  are  within  the  boundaries  of  empirical  shear  zone 
ratio.  Figures  (7.31)  and  (7.38)  show  the  velocity  profiles  within  the  shear  zone  in  the 
sliding  direction  for  slider  bearings.  It  is  expected  that  the  pressure  gradient  in  the  exit 
region  would  be  large  due  to  the  small  dimensions  of  the  shear  zone  ratio.  Figures  (7.32), 
(7.33),  and  (7.39)  show  the  film  temperature  distribution  for  the  slider  bearing.  Unlike 
the  results  for  the  journal  bearing  they  show  some  convection  effect  in  both  layers.  The 
high  temperature  lines  around  exit  zone  in  the  shear  zone  is  due  to  the  thin  shear  zone  in 
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that  region.  Heat  generation  takes  place  in  the  shear  zone  which  is  located  in  the  center 
layer  of  the  film  and  heat  is  carried  away  by  convection  as  well  as  conduction  (not 
considered  here  as  the  main  cause  of  the  temperature  reduction  on  the  solid  surface). 

A parametric  investigation  of  the  effect  of  film  thickness  on  the  shear  zone  was 
undertaken  for  the  case  of  a parallel  slider  with  B=L=3.0  inch  using  SAE  10  oil  at  100  ”F 
inlet  temperature  and  speed  U=10  fi/sec.  Film  thickness  values  from  0.002  to  0.1  inch 
were  considered.  The  shear  yield  stress  parameter  used  in  oil  cases  was  assumed  to  be 
fo  = 15.0  . The  maximum  pressure,  the  average  pressure  and  hs/h  at  the  inlet,  center  and 


exit  region  are  given  in  Table  11.  The  Sommerfelt  number  (S 
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BU)  are  calculated  for  each  case  and  the  empirical  h/h  was 


also  obtained  from  Equation  (7.58).  It  can  be  seen  from  the  table  that  the  shear  zone 
thickness  from  the  proposed  analysis  reaches  can  asymptotic  value  of  0.065942  which  is 
close  to  the  value  0.063662  calculated  from  the  empirical  equation.  It  can  also  be  seen 
that  the  proposed  method  shows  that  h/h  at  the  center  and  exit  region  increases  with 
increasing  h which  is  consistent  with  the  results  of  the  empirical  equation.  Figure  (7.40) 
shows  a plot  the  shear  zone  for  the  considered  case  as  obtained  from  the  proposed 


approach. 


7.9  Summary  and  Conclusion 

The  concept  of  the  shear  zone  as  determined  by  the  shear  yield  stress  parameter  is 
introduced  to  explain  the  effect  of  the  rheological  properties  of  the  lubricant  on  thin  film 
behavior  in  tribology.  The  shear  zone  is  hypothesized  to  be  composed  of  two  layers 
where  anti-symmetric  velocity  behavior  occurs  in  the  central  region  of  the  lubricating 
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film.  Interfacial  boundaries  between  the  shearing  fluid  and  the  adjacent  solid  like  fluid 
are  determined  from  the  points  where  the  resultant  shear  stress  is  less  than  yield  value. 
This  approach  produces  thin  shear  zone  at  the  central  region  of  the  film  where  the 
lubricant  undergoes  Newtonian  behavior  at  two  layers  with  mirror  image  dimensions 
about  the  central  plane  with  anti-symmetric  velocity  behavior.  The  parametric 
investigation  performed  using  different  rheological  properties  can  be  applied  to  any  thin 
film  lubrication  analysis  in  the  laminar  regime  but  does  not  appear  to  b applicable  to 
conditions  of  very  slow  running  speed  (N<1000  rpm). 


Table  1 1 : A Parametric  Investigation  of  the  Effect  of  Film  Thickness  on  the  Shear  Zone 
and  Summarized  Results  Based  on  the  Empirical  Formula  (7.58)  for  Evaluating  the  Shear 
Zone  Thickness. 


1 h (in) 

(h,/h)cal 

Pmax(psi) 

Pave(psi) 

s 

Pe(lO^) 

(hs/h)emp 

Inlet 

Center 

Exit 

0.002 

0.110924 

0.053987 

0.002581 

279.8932 

123.4971 

0.170199 

3.14496 

0.033312 

0.004 

0.083831 

0.056191 

0.032094 

33.1293 

13.8187 

0.380265 

3.14496 

0.052911 

0.006 

0.074351 

0.058484 

0.046004 

6.69668 

2.7692 

0.843368 

3.14496 

0.062571 

0.01 

0.068634 

0.060871 

0.056191 

0.66574 

0.2759 

3.047348 

3.14496 

0.063662 

6.02 

0.065942 

0.063356 

0.063356 

0.02341 

0.009989 

2 1 .04223 

3.14496 

0.063662 

6.1 

0.065942 

0.065942 

0.065942 

0.0000 1 

4.464  e-6 

1883.641 

3.14496 

0.063662 

Note:  B=L=3.0  inch,  SAE  10,  U=10  fit/sec,  Time=0.1  sec,  fa  = 15.0 
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Figure  7. 1 Shear  stress  vs.  Shearing  rate  for  Bingham  body;  Xq  is  the  yield  stress, 

Pseudo-plastic  or  shear-thinning,  and  Dilatant  or  shear-thickening.  Dashed  lines 
show  Newtonian  behavior  [79]. 
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Figure  7.2  Rigid  layers  attached  to  upper  and  lower  surfaces  for 
explaining  the  Non-Newtonian  behavior  in  fluid  film 
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Figure  7.3  (a)  Finite  lubrication  system  configuration  and  geometry 

(b)  The  computational  domain 
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Figue  7.4  Program  flow  chart  of  THD  analysis  using  the  Non-Newtonian 
velocity  expression  of  the  Rhelogical  behaviors 


Maximum  Pressure,  psi 
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Figure  7.5  Journal  bearing  pressure-moving  speed  characteristics 
D=3.154  in,  L=2  in,  C=0.0134  in,  SAE  10,  Tj„=128  °F,  8=0.9 

N=1050  rpm  : Transient  time  0.3  sec, 

N=2100  rpm  : Transient  time  0.5  sec,  Tq=15 


Maximum  Pressure,  psi 
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Figure  7.6  Journal  bearing  pressure-moving  speed  characteristics 
D=3.154  in,  L=2  in,  C=0.0134  in,  SAE  30,  Tj„=107  °F,  s=0.6 

N=1000  rpm  : Transient  time  0.1  sec,  Tq=8 
N=2000  rpm  : Transient  time  0.05  sec,  Xq=8 


Maximum  Pressure,  psi 
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Figure  7.7  Journal  bearing  pressure  history  for  various  yield 
shear  stress  parameter 


D=3.154  in,  L=2.0  in,  C=0.0134  in,  SAE  10,  Tj„=128  °F,  N=2100  rpm. 


8=0.9 


Maximum  Pressure,  psi 
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Figure  7.8  Journal  Bearing  pressure  history  for  various  shear  yield 
stress  parameters 

D=3.154  in,  L=2.0  in,  C=0.0134  in,  SAE  50,  Tj^=160  ®F,  N=2400  rpm, 
8=0.87 


Yield  shear  stress  parameter 
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Figure  7.9  Shear  zone  ratio  (hg/h),  based  on  the  THD  analysis 
with  various  shear  yield  stress  parameters. 
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Figure  7.10  Shear  zone  ratio  along  film  flow  direction 
D=3.154  in,  L=2.0  in,  C=0.0134  in,  s=0.9,  Tj„=128  ®F, 
SAE  10,  Time=0.5  sec,  N=2100  rpm 


Dimensionless  film  thickness 
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Figure  7.1 1 Shear  zone  ratio  along  film  flow  direction 
D=3.154  in,  L=2.0  in,  C=0.0134  in,  s=0.87,  Tjj^=160  °F, 
SAE  50,  Time=0.2  sec,  N=2400  rpm 


Dimensionlessshear  zone  location 
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Figure  7.12  Shear  zone  ratio  along  film  flow  direction 
D=3.154  in,  L-2.0  in,  C=0.0134  in,  e=0.6,  Tj^=107  ®F, 
Time=0.1  sec,  N=2000  rpm,  SAE  30 


Dimensionless  shear  zone,  h /h(x/B) 
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Figure  7.13  Velocity  profile  in  the  shear  zone  along  the 
mid  plane  of  the  journal  bearing  in  the  direction  of  motion 
D=3.154  in,  L=2.0  in,  C-0.0134  in,  s=0.87,  SAE  50, 

Tjn=160  °F,  N=2400  rpm.  Transient  time=0.2  sec,  T^=15 
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Figure  7.14  Pressure  distribution  for  various  shear  yield  parameters 
along  the  center-line  of  the  journal  bearing  in  the  direction  of  sliding 
motion,  D=3.154  in,  L=2.0  in,  C=0.0134  in,  8=0.6,  N=2000  rpm, 

SAE  30,  Tjjj=107  °F,  time=0.1  sec 


Maximum  pressure,  psi 
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Figure  7.15  Pressure  distribution  for  various  bearing  running  speed  along 
the  center  line  of  the  journal  bearing  in  the  direction  of  sliding  motion 

D=3.154  in,  L=2.0  in,  C=0.0134  in,  8=0.9,  SAE  10,  Tjj^=128  °F 


Moving  layer  Shear  zone  Stationary  layer 
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Figure  7.16  Dimensionless  film  temperature  distribution  at  the 
journal  bearing  mid  plane  in  the  sliding  direction 
D=3.154  in,  L=2.0  in,  C=0.0134  in,  8=0.9,  SAE  10 

Tjjj=128  °F,  Time=0.5  sec,  Tq=15,  N=2100  rpm 
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Figure  7.17  Dimensionless  film  temperature  distribution  at  the 
journal  bearing  mid  plane  in  the  sliding  direction 

D=3.154  in,  L=2.0  in,  C=0.0134  in,  8=0.6,  SAE  30,  Tij^=107°F, 

Time=0.05  sec,  Tq=7,  N=2000  rpm 


Moving  layer  Shear  zone  Stationary  layer 
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Figure  7.18  Dimensionless  film  temperature  distribution  at  the 
journal  bearing  mid  plane  in  the  sliding  direction 


D=3.154  in,  L=2.0  in,  C=0.0134  in,  e=0.6,  SAE  30,  Tj^=107  °F, 
Time=0.1  sec,  Tq=8,  N=1000  rpm 
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Figure  7.19  Slider  bearing  pressure-moving  speed  characteristics 
B=L=3.0  in,  h^- =0.004  inch,  h^,,=0.008  inch,  SAE  10,  T- =100  °F, 

’ mm  ’ max  ’ ’ m ’ 

time=0.02  sec,  Tq=10.0 
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Figure  7.20  Maximum  pressure  and  temperature  history  for 
various  runnung  speeds,  B=L=3.0  in,  =0.004  in,  h =0.008  in, 

To=10.0,  SAE  10,  Tj^=100 
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Figure  7.21  Slider  bearing  maximim  pressure  and  temperature 
history  for  various  shear  yield  stress  parameters, 

B=L=3.0  in,  h . =0.004  in,  h =0.008  in,  U=40  ft/sec, 

SAE  10,  T.  =100  ®F,  Time=0.2  sec 

’ in  ’ 


Pressure,  psi 
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Figure  7.22  Pressure  distribution  along  the  center  line  of  the 
slider  bearing  in  the  directionof  sliding. 


SAE  10,  h^.  =0.004  in,  h^,,=0.008  in,  T- =100  ®F, 

’ mm  ’ max  ’ in  ’ 

Time=0.02  sec,  U=40  fl/sec,  B=L=3.0  in 
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Figure  7.23  Pressure  distribution  along  the  center  line  of  the 
slider  system  in  the  direction  of  sliding, 


o, 


SAE  10,  h^i„-0.004  in,  h^^^=0.008  in,  Tj„=100  "F, 


U=20  ft/sec,  B=L=3.0  in. 
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Figure  7.24  Shear  zone  ratio  with  two  other  layers  for  various 
shear  yield  stress  parameters,  B=L=3.0  in,  hjjjjj^=0.004  in, 

h =0.008  in,  Tin=l  00  "F,  U=40  ft/sec,  and  SAE  1 0 

max  ’ ’ ’ 
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Figure  7.25  Shear  zone  ratio  (hj/h),  based  on  thermohydrodynamic 

analysis  with  various  shear  yield  stress  parameters 
hj^jjj=0.004  in,  hj^^j^^O.OOS  in,  B=L=3.0  in, 

Tjn=100  °F,  U=40  ft/sec,  SAE  10 


Dimensionless  film  thickness 


Figure  7.26  Shear  zone  ratio  along  the  film  flow  direction 
B=L=3.0  in,  h„,.  =0.004  in,.  h„,,  =0.008  in,  T- =100  °F, 

’ mm  ’ max  ’ in  ’ 

Time=0.2  sec,  U=40  ft/sec,  SAE  10,  Xo=10 
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Figure  7.27  Shear  zone  ratio  (hg/h),  based  on  thermohydrodynamic 

analysis  with  various  shear  yield  stress  parameters 
h«iin“0.004  in,  h,„„  =0.004  in,  B=L=3.0  in, 

mm  ’ max  ’ ’ 

Tin=100  °F,  U=20  fVsec,  SAE  10 


Dimensionless  shear  zone  thickness 


Figure  7.28  Shear  zone  ratio  along  film  flow  direetion 
B=L=3.0  in,  h„. -0.004  in,.  -0.008  in,  T- -100  °F, 

’ mm  ’ max  ’ in  ’ 

Time-0.2  sec,  U-20  ft/sec,  SAE  10,  Tq=10 
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Figure  7.29  Shear  zone  ratio  (hg/h),  based  on  thermohydrodynamic 

anaysis  for  various  shear  yield  stress  parameters 
h^- =0.004  in,  h^„^=0.008  in,  B=L=3.0  in, 

min  ’ max  ’ ’ 

Tj„=100  °F,  U=10  ft/sec,  Time=0.1  sec,  SAE  30 


Dimensionless  shear  zone 


Figure  7.30  Shear  zone  ratio  along  the  film  flow  direction 
B=L=3.0  in,  h^i„=0.004  in,.  h^^=0.008  in,  Tj„=100  °F, 
Time=0.1  sec,  U=10  ft/sec,  SAE  30,  Xo=10 


Shear  Zone,  h /h(x/B) 
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Figure  7.31  Velocity  profile  in  the  shear  zone,  h , along  the 

mid-plane  of  the  slider  bearing  in  the  direction  of  motion 
B=L=3.0  in,  li- =0.004  in,  h^,  =0.008  in,U=40  ft/sec, 

’ mm  ’ max  ’ ’ 

SAE  10,  Xq=10.0,  Tjj^=100  ”f,  time=0.02  sec 


Moving  layer  Shear  zone  Stationary  layer 
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Figure  7.32  Dimensionless  film  temperature  distribution  along  the 
central  line  of  the  slider  bearing  in  the  direction  of  film  moving, 

B=L=3.0  in,  h^j„=0.004  in,  10,  Tjjj^lOO  °F, 

t'q=10.0,  U=40  ft/sec,  time=0.02  sec 


Moving  layer  Shear  zone  Stationary  layer 

(dimensionless)  (dimensionless)  (dimensionless) 
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Figure  7.33  Dimensionless  film  temperature  distribution  along  the 
central  line  of  the  slider  bearing  in  the  direction  of  film  moving, 


o. 


B=L=3.0  in,  li.  =0.004  in,  li,  =0.008  in,  SAE  30,  T- =100  F, 
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Figure  7-34  Parallel  slider  bearing  maximum  pressure  and 
temperature  history  for  various  parallel  clearances 
B=L=3.0  in,  h=0.004  in  (Tq=13),  0.006  in  (Xq=15),  and  0.008  in, 

(T(j=27),  U=10  ft/sec,  SAE  30,  Tj^=100  "F,  Time=0.1  sec 
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Figure  7.35  Pressure  distribution  along  the  central  line  of  the 
parallel  surface  bearing  in  the  direction  motion, 

B=L=3.0  in,  SAE  30,  Tj„=100  ®F,  U=10  ft/sec,  h=0.006  in, 
Time=0.1  sec 


Pressure,  psi 
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Figure  7.36  Pressure  distribution  along  the  central  line  of 
the  slider  bearings  in  the  direction  of  motion,  B=L=3.0  in, 

SAE  30,  Tj^=100  °F,  U=10  ft/sec.  Transient  time=0.1  sec 

Case  J - Parallel  slider,  h=0.004  in,  Tq=13 
Case  K - Parallel  slider,  h=0.006  in,  fQ=15 
Case  L - Parallel  slider,  h=0.008  in,  Xq=27 
Case  I - Tapered  slider, 
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Figure  131  Shear  zone  ratio  (h„/h),  based  on  thermohydrodynamic 
analysis  for  various  shear  yield  stress  parameters,  {=x  hQ/(p  U)}, 

B=L=3.0  in,  SAE  30,  Tj^=100  °F,  U=10  ft/sec,  h=0.006  in,  Time=0.1  sec 
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Figure  7.38  Velocity  profile  in  the  shear  zone,  h„/h,  along  the 

mid-plane  of  the  parallel  slider  bearing  in  the  direction  of  motion 
B=L=3.0  inch,  hjjjjj^=0.006  inch,  U==10  ft/sec,  SAE  30,  Tq=15.0, 

Tin=100  °F,  time=0.1  sec 
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Figure  7.39  Dimensionless  film  temperature  distribution  along  the 
central  line  of  the  parallel  slider  bearing  in  the  direction  of  film  moving 

B=L=3.0  in,  h=0.006  in,  SAE  30,  Tjj,=100  ®F,  Xo=15.0, 

U=10  ft/sec,  Time=0.1  sec 


h„/h.  Thickness  of  shear  zone  ratio 
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Figure  7.40  The  trend  of  the  shear  zone  ratio  with  various  film 
thickness  along  the  film  flow  direction  in  the  parallel  slider 


B=L=3.0  in,  Tjj,=100  °F,  Time=0.1  sec,  U=10  ft/sec. 


SAE  10,fo=15 


CHAPTER  8 

SUMMARY,  CONCLUSIONS  AND  RECOMMENDATIONS 


8.1  Summary  and  Conclusions 

The  main  goal  of  the  thesis  is  to  develop  a mathematical  procedure  for  more 
accurate  prediction  of  the  film  thickness  as  well  as  the  velocity,  shear  and  pressure 
distribution  in  fixed  geometry  tribological  systems.  In  response  to  this  goal  numerical 
solutions  are  presented  for  flow  condition  on  a free  surface,  free  surface  with  constraint 
and  thermohydrodynamic  lubrication  between  two  surfaces.  The  effect  of  surface 
tension,  speed,  system  geometry,  and  lubricant  rheological  properties  on  the  velocity, 
temperature  and  pressure  distribution  in  the  film  are  investigated.  The  following 
summarizes  the  findings  of  the  reported  study: 

Some  general  analysis  for  the  free  surface  generated  by  the  withdrawal  process  is 
presented  as  an  extension  of  the  Landau-Levich  theory.  The  film  thickness  of  the 
meniscus  is  found  to  affected  primarily  by  Capillary  number,  (Ca),  at  low  speeds  and  also 
by  Reynolds  number,  (Re),  at  high  speeds.  In  the  domain  of,  Ca  < 10,  the  results  for 
Newtonian  liquids  on  inclined  plates  are  obtained.  Beyond  Ca  = 10,  it  is  important  to 
consider  the  effect  of  the  inertia  force  in  the  analysis  of  the  free  surface  by  the  withdrawal 
process.  The  analytical  solution  was  obtained  as  an  extension  of  the  flat  plate  case  by 
applying  geometric  analogy  to  a rotating  disk.  The  depth  of  the  submerged  portion  of  the 
disk  in  the  oil  container  was  found  to  have  significant  effect  on  the  film  thickness 
generated  by  surface  tension.  The  numerical  solutions  for  predicting  the  withdrawal  film 
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thickness  on  the  rotating  disk  has  been  shown  to  be  stable  for  Ca  <10  and  in  good 
agreement  with  experiments  using  a 6 inch  diameter  disc  rotating  at  speeds  below  180 
rpm.  The  film  thickness  and  mass  flow  rate  on  the  rotating  disk  can  consequentially  be 
used  to  determine  the  boundary  conditions  in  the  lubrication  problem. 

A general  theoretical  treatment  of  steady,  two-dimensional,  Newtonian  film 
splitting  flows  between  a rotating  disc  and  a stationary  surface  has  been  presented.  The 
simple  mathematical  formulation  based  on  the  lubrication  model  can  predict  the 
dimensionless  film  thickness  to  be  formed  after  the  outlet  zone.  It  is  found  to  be  a 
function  of  the  two  boundary  points  and  the  gravity  number.  To  determine  the  two 
boundary  points,  one  of  the  models  used  considers  the  Swift-Stieber  (Reynolds)  boundary 
condition  at  the  exit  region  where  the  lubrication  flow  splits.  The  other  model  considers 
the  continuity  boundary  condition  governed  by  the  geometrical  parameter  Hq/R 
(minimum  gap/disc  radius),  rotating  speed,  and  Ca  using  the  mathematical  techniques 
developed  by  Coyne  and  Elrod.  The  splitting  film  thickness  is  expressed  by  the 
dimensionless  flow  rate  and  the  splitting  position.  This  position  is  determined  from  the 
entrance  position  which  is  based  on  the  previous  results  of  the  withdrawal  process.  The 
major  point  is  that  the  Swift-Stieber  (Reynolds)  boundary  condition  is  not  suitable  for 
dealing  with  the  existence  of  a negative  pressure  before  the  separation  and  the  effect  of 
surface  tension  which  is  significant  at  low  speeds.  At  a critical  speed,  there  exists  an 
interaction  between  the  pressure  distribution  before  the  meniscus  region  and  the 
interfacial  pressure  controlled  by  surface  tension.  Beyond  that  speed  the  phenomenon 
can  be  simply  modeled  by  Swift-Stieber  (Reynolds)  boundary  condition  which  is  used  in 
the  analysis  of  slider  and  journal  bearing  systems. 
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A numerical  method  is  developed  for  the  THD  lubrication  of  slider  and  journal 
bearings  considering  both  the  fiill  shear  zone  and  the  partial  shear  zone  taking  into 
consideration  the  thermal  effect.  The  computational  technique  requires  the  simultaneous 
solution  of  the  Newtonian  hydrodynamic,  the  viscosity  variation  and  the  energy 
equations.  The  results  include  the  lubricant  pressure  profile,  film  thickness,  and 
temperature  distribution  within  the  lubrication  zone.  The  analysis  leads  to  the  following 
main  conclusions: 

1.  The  results  of  the  Newtonian  behavior  in  the  full  shear  zone  show  that  the 
temperature  rise  in  the  lubrication  zone  is  proportional  to  the  sliding  speed,  but  the 
thermal  effects  are  not  significant  enough  to  reduce  the  pressure  in  the  lubrication  region. 

2.  The  major  component  of  heat  generation  under  sliding  conditions  is  by  shear  and 
for  cooling  is  by  conduction. 

3.  The  relationship  between  the  operating  speed  and  film  pressure  in  a fixed  geometry 
bearing  as  obtained  from  both  the  Newtonian  and  Rheological  models  is  found  to  follow 
the  square  root  dependency  of  pressure  on  speed  as  proposed  by  Seireg. 

4.  The  use  of  the  Bingham  model  in  fluid  film  lubrication  is  introduced.  The  non- 
Newtonian  fluid  under  consideration  differs  from  the  Newtonian  fluids  in  the  existence  of 
a narrow  central  layer  in  the  lubrication  zone  where  the  velocity  gradients  are  exceedingly 
high.  From  this  consideration  an  anti-symmetric  velocity  model  for  the  Newtonian  fluids 
was  presented. 

5.  The  Bingham  model  with  anti-symmetric  velocity  profiles  in  the  shear  zone  was 
incorporated  with  the  Reynolds  equation  and  the  energy  equation  for  THD  lubrication 
analysis. 

6.  The  shear  zone  was  determined  from  the  condition  where  the  shear  stress  exceeds 
the  shear  yield  stress  parameter  which  defines  the  interfacial  boundary  between  the  solid 
and  fluid  behavior.  The  computational  method  was  used  to  evaluate  the  shear  yield  stress 
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parameter,  f o , which  gives  values  for  the  film  the  maximum  pressure  corresponding  to 
the  experimental  data. 

7.  The  results,  show  that  the  present  model  can  be  applied  for  lubrication  conditions 
with  running  speeds,  U > 10  ft/sec  in  the  slider  bearings  and  in  the  journal  bearings 
corresponds  to  N>1000  rpm.  It  is  shown  that  the  shear  yield  stress  parameter  as 
determined  for  a journal  bearing  or  slider  bearing  is  applicable  for  the  same  geometric 
parameters  at  different  running  speed. 

8.  Except  very  low  speeds  (U  < 10  ft/sec  or  N < 1000  rpm)  and  very  small  clearance 
hmin  < 0.001  inch)  the  range  for  dimensionless  shear  yield  stress  parameter  in  investigated 
bearings  is  7 ~ 15. 

9.  The  calculated  thickness  of  the  shear  zone  along  the  lubrication  flow  direction  was 
compared  with  the  empirical  relationship  [28]  for  calculating  shear  zone  dimensions  as  a 
function  of  the  product  of  Sommerfeld  number  (S)  and  Peclet  number  (Pe).  The 
calculated  shear  zone  dimensions  by  the  shear  yield  stress  parameter  are  of  the  same 
order  of  magnitude  as  those  calculated  from  the  empirical  relationship. 

10.  The  results  of  the  parametric  study  for  the  effect  of  the  film  thickness  in  the  shear 
zone  in  the  parallel  slider  bearing  show  the  same  trend  and  approximately  the  same 
asymptotic  value  for  h^/h  as  predicted  by  the  empirical  equation. 

8.2  Recommendations 


1.  Needless  to  say  further  and  more  detailed  studies  need  to  be  performed  for  full 
understanding  of  the  withdrawal  problems  including  their  applications  in  lubrication. 
More  attention  should  be  focused  on  the  effect  of  surface  tension  on  the  film  thickness  in 
the  free  surface  and  on  the  film  flow  behavior  in  the  clearance  between  the  rotating  disc 
and  stationary  surface.  The  effect  of  surface  tension  acting  in  the  entrance  region  where 
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lubrication  begins  which  is  known  to  be  influenced  by  the  rheological  behavior  of 
lubricants  and  the  modified  Bingham  model  in  the  lubrication  region  should  be  extended 
by  considering  the  effect  of  the  surface  tension. 

2.  The  present  research  has  provided  further  insight  for  better  understanding  of  THD 
lubrication  under  the  considered  conditions.  There  are  many  opportunities  to  continue 
research  in  each  of  the  problem  areas  investigated.  The  derived  simple  form  of  flow 
under  THD  lubrication  can  be  readily  extended  to  incorporate  more  Non-Newtonian 
effects. 

3.  Broader  experimental  confirmation  for  the  theoretical  models  is  also  necessary. 
The  results  obtained  from  two  experimental  slider  bearing  and  two  journal  bearing  cases 
should  be  checked  for  a wider  range  of  bearing  geometries  and  speeds. 

4.  A number  of  problems  of  thermohydrodynamic  lubrication  under  special 
conditions  have  been  studied.  These  include  cases  where  the  thermal  effect  is  important. 
The  Reynolds  equation,  energy  equation,  and  a semi-empirical  equation  have  been 
coupled  to  solve  the  THD  problem.  However,  in  addition  to  the  rhelogical  phenomenon 
with  an  adiabatic  high  shear  rate  band  in  the  lubrication  region,  heat  conduction  and 
cavitation  effects  must  be  also  considered  in  further  studies. 

5.  The  shear  yield  stress  is  shown  in  this  study  to  have  a major  influence  on  the 
shear  zone  and  the  corresponding  behavior  of  the  lubricating  film.  Developing  reliable 
data  on  its  magnitude  for  different  oils  at  different  temperatures,  pressure  and  shear  rates 
would  be  necessary  for  the  prediction  of  film  behavior. 

6.  When  reliable  shear  yield  stress  data  are  available,  it  would  be  interesting  to  study 
the  effect  of  film  thickness  and  speed  on  the  maximum  pressure  using  both  the  full  film 
shear  and  partial  shear  zone  assumptions. 
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